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ABSTRACT 


The  second  AFSWC  Hydrodynamic  Conference  was  held  at  Kirtland  Air 
Force  Base  for  the  purpose  of  discussing  numerical  techniques  which  have 
become  an  important  and  widely  used  tool  for  solving  fluid  flow  problems. 

A  large  number  of  the  papers  presented  dealt  with  the  problems  of  the  finite 
difference  analogs  of  the  differential  equations  of  motion.  A  few  papers 
discussed  work  toward  analytic  solutions  of  these  equations.  Topics  in¬ 
cluded  hydrodynamics,  magnetohydrodynamics,  radiation  transport,  and 
solid  material  motion.  (U) 

The  Proceedings  are  published  in  two  parts;  Part  1  is  unclassified  and 
Part  II  is  classified.  (U) 
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WELCOMING  ADDRESS 

Col  David  R,  Jones 
Chief,  Physics  Division,  AFSWC 


I  would  like  to  welcome  you  to  Kirtland  on  behalf  of  General  McCorkle, 
the  Research  Directorate,  and  the  Physics  Division. 

This  conference  can  very  well  be  considered  as  the  second  AFSWC 
Hydrodynamics  Conference.  The  first  was  held  in  March  of  1960.  At  that 
time  we  were  primarily  concerned  with  the  X-ray  effects  of  a  nuclear  burst, 
and  needed  to  arrive  at  new  levels  of  common  understanding  before  embarking 
on  an  integrated  large-scale  in-house  and  contractual  research  program. 

This  year  we  are  not  concerning  ourselves  with  an  individual  program, 
but  are  attempting  to  discuss  one  aspect  of  a  number  of  the  Physics  Division's 
programs.  This  aspect  is  hydrodynamics.  We  have-chosen  to  concentrate  on 
numerical  methods  and  techniques  of  hydrodynamics,  but  have  included  some 
topics  dealing  with  the  physics  involved  and  w-ith  non-numerical  methods  of 
solution  to  these  problems. 

The  overall  purpose  of  this  meeting  should  be  toward  the  exchange  of 
methods  and  techniques  of  approaching  and  solving  various  problems.  We 
hope  that  the  idea  of  a  working  meeting  can  be  kept  in  mind  and  that  impromptu 
discussion  of  unresolved  topics  will  flow  freely. 

Again,  I  repeat  my  welco.me  to  Kirtland,  and  offer  our  services  to  you. 
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ON  THE  TREATMENT  OF  RAREFACTION 
USING  A  DISSIPATIVE  HYDRODYNAMICS  CODE 

Gunning  Butler.  Jr.  ,  and  Daniel  M.  Young 
Boeing  Airplane  Company 


One  of  the  methods  being  used  to  study  the  response  of  materials  to  the 
X-ray  impulse  is  slapping  the  material  with  a  thin  flying  plate.  The  Boeing 
group  working  on  this  problem  is  accelerating  the  flyer  plates  by  discharging 
a  condenser  bank  through  a  thin  metallic  foil  and  using  the  pressure  de¬ 
veloped  in  the  exploded  foil  to  accelerate  the  flyer.  After  the  plate  has  moved 
a  sufficient  distance  it  is  essentially  at  initial  density  and  moving  with  con¬ 
stant  velocity.  When  the  plate  impinges  on  the  target  material,  a  square 
pressure  pulse  is  generated  in  the  material  with  a  width  equal  to  the  time  for 
the  resulting  shock  to  travel  back  through  the  flyer  plate  and  the  rarefaction 
from  the  free  surface  to  travel  back  to  the  contact  surface.  As  the  shock 
progresses  into  the  target  material  it  will  develop  into  a  normal  triangular 
shock  because  of  rarefactions  from  the  surface  relieving  the  back  of  the  pulse 
and  degrading  the  peak  pressure. 

In  order  to  understand  the  details  of  these  processes  and  the  subsequent 
details  of  the  flow,  a  hydrodynamic  code  was  developed  at  Boeing.  The  code 
is  based  on  the  Richtmyer- Von  Neumann  method  and  the  details  of  the  code 
are  reported  in  the  proceedings  of  the  Air  Force  Special  Weapons  Center 
Hydrodynamics  Conference  (AFSWC  TR-60-1Z)  which  was  held  here  about  a 
hear  ago.  ^ 

The  manner  in  which  a  dissipative  hydrodynamics  code  such  as  that  of 
Von-Neumann  and  Richtmyer^  or  Ludford,  Polachek .  and  Seeger^  handles 
one-dimensional  shocks  is  well  established  and  fully  reported  in  the  litera¬ 
ture.  The  usual  set  of  basic  difference  equation  is 
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For  reference,  the  pertinent  phenomenological  features  are  shown  in  figure  1, 
which  is  a  synoptic  view  of  a  shock  in  terms  of  fluid  particle  velocity  as  cal¬ 
culated  by  the  Boeing  hydrodynamics  code  and  compared  with  the  idealized 
shock  in  a  perfect  fluid.  The  independent  variable  here  is  the  Lagrangian 
coordinate.  One  notes  the  spreading  of  the  shock  front  by  the  code  and  the 
oscillations  or  ripples  induced  in  the  flow  after  passage  of  the  shock.  Figure 
2  shows  the  pressure  time  history  of  a  mass-zone  after  the  shock  has  passed. 
It  is  well  known  that  as  the  dissipation  terms  in  the  equations  are  increased  in 
magnitude,  the  amplitude  of  the  oscillations  is  decreased  and  the  shock  thick¬ 
ness  is  increased. 

Originally  it  was  decided  to  start  the  calculations  by  dumping  the  energy 
based  on  the  measurements  of  the  current  and  voltage  versus  time  into  the 
foil.  However,  it  was  not  possible  to  treat  this  situation  with  sufficient 
accuracy.  The  results  of  dumping  the  energy  into  the  foil  are  extremely 
sensitive  to  the  equation  of  state  of  the  foil.  A  series  of  runs  was  done  in 
which  the  foil  was  treated  as  an  ideal  gas  and  the  peak  pressures  and  peak 
free-surface  velocities  varied  quite  widely  with  a  small  change  in  gamma. 

This  effort  was  then  dropped  because  of  lack  of  knowledge  of  the  foil  equation 
of  state.  At  the  present  time  the  input  is  taken  as  the  flyer  plate  at  normal 
density  and  constant  velocity  at  the  time  it  impinges  on  the  target,  because 
measurements  of  the  flyer  plate  velocity  are  readily  measurable  and  quite 
reproducible.  Also  by  using  the  flyer  plate  in  motion  as  input,  the  equation 
of  state  of  the  foil  need  not  be  known  and  only  the  equation  of  state  of  the  flyer 
over  a  relatively  small  pressure  range  is  required. 

The  Boeing  group  is  particularly  interested  in  the  mechanism  of  spall. 
Spall  occurs  ,vhen  a  material  is  in  tension.  For  this  reason  it  is  important 
to  know  how  a  Von-Neumann-type  code  treats  rarefactions.  If  these  codes  do 
not  handle  rarefactions  properly  then  there  are  problems  in  determining  the 
location  of  the  maximum  tension  and  its  absolute  magnitude.  For  this  reason 
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a  series  of  numerical  experiments  on  centered  rarefactions  has  been  conducted 
at  Boeing.  Mr.  Daniel  M.  Young  will  now  present  the  results  of  these  studies. 

At  the  suggestion  of  Dr.  John  G.  Trulio  of  the  Lawrence  Radiation  Lab¬ 
oratory,  Livermore,  the  Boeing  group  has  undertaken  what  can  best  be  de¬ 
scribed  as  a  series  of  numerical  experiments  on  simple  centered  waves  using 
the  Boeing  hydrodynamics  code.  The  treatment  of  these  types  of  flow  by  such 
codes  is  of  some  interest,  but  does  not  seem  to  have  been  reported  in  the 
literature,  at  least  in  the  available  recent  and  near-recent  sources.  For 
normal  materials,  simple  centered  waves  are  rarefactions  and  are  repre¬ 
sented  by  similarity  solutions  of  the  equations  of  motion;  i.  e.  ,  the  equations 
are  reduced  to  a  set  of  ordinary  differential  equations  which  can  be  integrated 
by  standard  methods.  For  certain  equations  of  state  such  solutions  are  ex- 
pressable  in  closed  form. 

Simple  centered  waves  are  of  interest  primarily  because  they  occur 
frequently  as  integral  parts  of  more  complicated  flow  fields  which  require 
a  hydrodynamic  code  for  solution.  The  manner  in  which  a  code  handles 
simple  centered  waves  in  cases  where  the  solution  is  known  should  be  in¬ 
dicative  to  some  extent  of  the  error  introduced  in  the  more  complicated 
flows  which  include  them.  This  is  particularly  important  when  one  is  trying 
to  look  at  spall  which  is  largely  determined  by  crossing  rarefactions.  The 
exact  similarity  solution  representing  a  simple  centered  wave  is  presented 
in  figure  3  which  is  a  synoptic  view  of  a  rarefaction  moving  into  a  medium  at 
rest.  The  material  in  this  case  has  a  Hooke's  law  equation  of  state.  Both 
pressure  and  particle  velocity  are  shown  as  functions  of  the  Lagrangian  co¬ 
ordinate. 

If  the  dissipation  terms  in  the  above  set  of  difference  equations  are  set 
to  zero  and  the  same  problem  is  calculated  with  the  hydrodynamics  code,  the 
results  are  those  shown  in  figure  4,  The  features  contrasting  the  analytical 
and  code  solutions  are  the  overshoot  in  the  flow  behind  the  wave  and  the 
spreading  of  the  wave  over  too  wide  a  region.  If  the  dissipation  terms  are 
included  in  the  calculation,  the  amplitude  of  the  overshoot  is  markedly  de¬ 
creased  but  the  spreading  is  increased,  as  is  indicated  in  figure  5  which 
shows  the  rarefaction  wave  both  with  and  without  dissipation.  As  in  the  case 
of  shock  compressions,  increasing  the  magnitude  of  the  dissipation  terms 
enhances  these  effects.  That  this  should  be  the  case  also  for  rarefactions  is 
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not  altogether  obvious,  since  for  such  flows  the  Von  Neumann- Richtmyer 
dissipation  is  negative.  The  inclusion  of  such  terms,  therefore,  corresponds 
to  a  physically  unreal  situation,  since  a  negative  dissipation  implies  violation 
of  the  second  law  of  thermodynamics.  It  should  be  recalled,  however,  that 
the  dissipation  terms  were  originally  introduced  into  the  equations  only  as  a 
mathematical  device,  not  as  an  attempt  to  simulate  physical  reality.  The 
magnitude  of  these  terms  required  to  allow  the  shock  calculation  to  proceed 
is  greatly  in  excess  of  the  real  physical  magnitude. 

Figure  6  is  a  wave  diagram  in  Lagrangian  coordinates  of  the  development 
of  a  simple  centered  wave  from  a  pressure  discontinuity  as  computed  by  the 
hydrodynamics  code  without  dissipation.  It  resembles  somewhat  grossly  the 
early  history  of  the  development  of  a  shock  compression  from,  a  pressure 
discontinuity,  shown  for  comparison  in  figure  7.  The  effect  of  the  dissipation 
term  is  shown  by  the  wave  diagrams  of  figure  8.  The  parameter  C  is  the 
constant  in  the  Richtmyer- Von  Neumann  dissipation  term.  Comparison  with 
the  wave  diagram  of  the  analytical  solution  indicates  that,  as  calculated  by 
the  code,  the  head  of  the  wave  initially  propagates  too  quickly  and  the  tail 
too  slowly.  Finally  the  wave  settles  down  to  a  stable  fan- shaped  structure 
with  all  parts  traveling  at  their  correct  speeds.  The  center  of  the  wave 
seems  to  travel  at  the  right  speed  from  the  beginning. 

There  seem,  then,  to  be  two  problems  associated  with  the  treatment  of 
simple  centered  waves  by  dissipative-type  hydrodynamic  codes:  (1)  there  is 
a  spurious  overshoot  and  subsequent  decaying  ripple  induced  on  '-.he  flow 
following  the  wave,  which  if  the  dissipative  terms  are  small  or  ignored 
completely,  can  be  of  significant  magnitude;  and  (2)  the  trajectories  of  the 
head  and  tail  of  the  wave  are  incorrect  because  they  are  initially  of  the  wrong 
slope,  a  situation  aggravated  by  increasing  the  dissipation. 

Inspection  of  the  details  of  the  difference  equations  at  the  contact  surface 
across  which  the  pressure  discontinuity  originally  exists  leads  one  to  suspect 
that  some  of  the  difficulty  is  the  result  of  the  way  in  which  the  velocity  of 
that  surface  is  computed.  The  usual  difference  equation  for  the  velocity  of 
a  mass  zone  boundary  is 
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where  the  nomenclature  is  essentially  that  of  Von  Neumann  and  Richtmyer. 
The  customary  way  of  handling  a  surface  is  to  use  this  same  form  of  differ¬ 
ence  equation.  Letting  P^^  denote  the  outside  pressure  on  the  boundary,  the 
equation  for  the  velocity  of  the  surface  is 
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where  the  boundary  has  been  assumed  to  be  on  the  right. 

Actually  the  right  hand  side  of  the  above  equation  (*)  represents  the 
secant  centered  one  quarter  of  a  mass  zone  from  the  surface.  For  simplicity 
consider  the  case  where  P^  is  zero.  Then  according  to  (*)  the  pressure  in 
the  first  zone  adjacent  to  the  surface  must  be  completely  relieved  before  the 
acceleration  of  the  surface  will  cease.  Thus  the  velocity  of  the  surface  will 
overshoot,  a  tension  will  develop  to  slow  it  down,  etc.  This  effect  will  then 
propagate  along  with  the  tail  of  the  wave.  The  early  history  of  the  surface 
velocity  is  shown  in  figure  It  is  interesting  to  note  that  when  dissipation 
is  included,  the  magnitude  of  the  initial  overshoot  is  decreased  as  it  should 
be  according  to  (*),  since  the  dissipation  term  will  always  be  opposite  in 
sign  to  the  pressure  difference  for  a  rarefaction  wave. 

When,  in  the  hydrodynamics  code,  the  rarefaction  is  initiated  by  the 
reflection  of  a  shock  wave  from  a  free  pressure  boundary,  it  is  not  a  simple 
centered  wave.  The  structure  of  the  shock  prevents  the  wave  from  being 
centered.  This  is  shown  in  figure  10,  which  is  a  wave  diagram  in  Eulerian 
coordinates  of  this  process.  However,  since  the  resulting  wave  is  not 
centered,  it  is  not  liable  to  such  violent  overshoots  in  its  formative  stages 
and  one  of  the  problems  has  been  alleviated.  The  cost  however,  is  a  larger 
error  in  the  trajectories  of  the  head  and  the  tail. 
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The  effect  of  such  errors  in  a  rather  simple  problem  is  illustrated  in 
figure  11.  The  problem  is  that  of  an  unstrained  flyer  plate  impacting  against 
a  semi-infinite  slab  of  identical  material.  The  quantities  of  interest  are  the 
distance  and  time  at  which  the  rarefaction  overtakes  the  shock,  and  the 
trajectory  and  shape  of  the  subsequently  decaying  shock  pulse.  The  sort  of 
error  that  is  prone  to  occur  is  too  short  a  catch-up  distance  after  which  the 
pulse  calculated  by  the  code  continues  to  fall  behind  because  of  its  early  loss 
of  driving  pressure.  The  effect  of  some  of  the  program  parameters  upon  the 
results  is  shown  in  figure  12.  The  parameter  to  which  the  catch-up  distance 
is  most  sensitive  is  the  dissipation  constant  C.  The  reason  is,  as  has  been 
shown,  that  the  smaller  the  dissipation  the  more  accurately  localized  are 
the  waves.  However,  simultaneously  the  amplitudes  tend  to  become  rather 
ill-defined  because  of  the  overshoots  and  subsequent  oscillations  in  the  flow 
regions  immediately  following  the  waves. 

It  seems,  then,  that  some  care  should  be  taken  in  interpreting  the  results 
of  problems  calculated  with  dissipative  hydrodynamics  codes  when  simple 
centered  waves  are  features  of  the  resulting  flows.  In  particular,  tension 
regions  following  rarefactions  should  be  examined  by  an  independent  method 
to  establish  their  validity.  Also,  if  possible,  catch-up  distances  should  be 
checked. 
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Figure  1.  Particle  velocity  vs.  zone  number 
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igure  4.  Comparison  of  analytical  and  code  solution  for  a  centered 
simple  wave  in  a  Hooke's  Law  material.  No  dissipation  in  code 
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Figure  6.  Early  history  of  development 
of  a  rarefaction  from  a  pressure  discontinuity 
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Figure  10.  Rarefaction  caused  by  shock  reflection 
at  a  free  pressure  boundary 
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- ANALYTICAL  SOLUTION 


Figure  11.  Degrading  of  shock  by  overtaking  rarefaction 
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ANALYTIC  METHODS  AND  APPROXIMATIONS  OF  MHD  PROBLEMS 

Dr.  J.  D.  Cole 

California  Institute  of  Technology 
and 

Dr.  C.  Greifinger 
The  RAND  Corporation 

I.  INTRODUCTION. 

In  recent  years,  largely  because  of  the  interest  in  controlled  thermo¬ 
nuclear  reactions,  many  devices  have  been  designed  to  accelerate  gases  to 
thermonuclear  energies  (~10  kev).  In  the  operation  of  most  of  these  devices, 
the  gas  is  driven  electromagnetically,  it  being  possible  in  this  manner  to 
achieve  much  higher  gas  velocities  than  by  mechanical  or  gas-dynamic  means. 

In  the  usual  mode  of  operation,  electrical  energy,  stored  in  a  charged 
capacitor,  is  delivered  to  the  apparatus  containing  the  gas.  The  gas  is 
rapidly  broken  down  by  the  applied  voltage,  and  currents  flow  through  it.  A 
magnetic  field  is  thereby  generated  which  interacts  with  the  currents  in  the 
gas  and  sets  the  gas  into  motion.  Generally,  the  currents  are  confined  to  a 
rather  thin  sheet  at  the  boundary  of  the  gas.  The  magnetic  field  then  acts  as 
a  piston,  pushing  the  boundary  inward  and  leaving  behind  a  region  of  vacuum. 
The  boundary  is  preceded  by  a  shock  wave  which  compresses,  ionizes,  and 
sets  into  motion  the  enveloped  gas. 

To  analyze  quantitatively  the  operation  of  such  devices,  it  is  obviously 
necessary  to  make  a  number  of  simplifying  approximations,  the  most  common 
of  which  are  the  following: 

(a)  idealization  of  the  geometry  of  the  apparatus, 

(b)  idealization  of  the  physical  properties  of  the  gas, 

(c)  neglect  of  the  reaction  of  the  gas  on  the  external  circuit,  and 

(d)  idealization  of  the  mechanical  properties  of  the  gas. 

The  geometrical  idealization  always  involves  the  neglect  of  end  effects 
and  occasionally  involves  other  idealizations,  one  of  which  will  be  illustrated 
below  in  the  analysis  of  the  inverse  pinch.  Such  approximations  do  not 
usually  impose  serious  limitations  on  the  validity  of  the  analysis,  but  merely 
limit  that  validity  to  some,  usually  large,  portion  of  the  apparatus. 
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The  idealization  of  the  physical  properties  of  the  gas  generally  consists 
in  the  assumptions  that  the  plasma  is  a  nonviscous,  non-heat-conducting, 
ideal  gas  and  that  it  has  infinite  electrical  conductivity.  These  approxima¬ 
tions  limit  the  validity  of  the  analysis  to  some  range  of  operating  conditions; 
fortunately,  under  normal  conditions,  the  devices  to  be  discussed  operate 
well  within  the  necessary  range.  This  point  is  best  illustrated  by  considering 
in  some  detail  the  assumption  of  infinite  conductivity. 

The  infinite  conductivity  approximation  is  equivalent  to  the  assumption 
that  the  currents  in  the  plasma  are  confined  to  an  infinitely  thin  sheet  on  the 
boundary  which  excludes  the  driving  magnetic  field  from  the  plasma. 

Actually  however,  because  of  the  finite  conductivity  of  the  plasma,  the  field 
diffuses  into  it.  The  diffusion  distance  &  at  a  time  t  is  roughly 


~  (-L_) 

p  0 


(1) 


where  a  is  the  conductivity  of  the  gas  and  |i  the  permeability  of  free  space. 
For  a  shock  of  roughly  constant  speed  c.  the  separation  A  of  the  shock 
and  the  current  carrying  contact  front  is 

A  ~  (c  -  u)t  (2) 

where  u  is  the  speed  of  the  contact  front.  The  density  ratio  across  the  shock 
is 
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where  the  last  equality  holds  for  a  very  strong  shock  in  a  perfect  gas  (y  1b 
the  usual  ratio  of  specific  heats).  For  a  real  gas  (with  partial  dissociation 
and  ionization),  (y  -  1)/(y  +  D  >  1/15,  so  that 
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The  approximation  should  be  valid  so  long  as  (taking  an  average  for  &  ) 
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If  this  is  to  apply  to  a  shock  traveling  a  distance  D,  than  t  ~  D/c  and 
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Thus,  a  magnetic  Reynolds  number  R^^  based  on  the  dimensions  of  the 
device  and  the  shock  speed  must  be  fairly  large, 

R-,  =  Li  dDc  >  100  ,  (7) 

M 

say.  As  a  typical  case,  for  a  Mach  20  shock  (10-ev  temperature)  progressing 

4  4 

through  cold  deuterium  c  ~  4  x  10  m/sec  and  d  ~  6  x  10  mho/m,  so  that 
the  inequality  required  by  equation  (7)  is  well  satisfied  for  dimensions  D 
greater  than  about  2  inches.  The  relevant  dimension  of  all  the  devices  to  be 
discussed  is  at  least  this  large,  so  that  the  approximation  may  be  expected 
to  apply  rather  well. 

The  approximation  of  neglecting  the  reaction  of  the  gas  on  the  external 
circuit  is  effected  by  replacing  the  external  circuit  by  a  boundary  condition. 
That  is,  the  current  through  the  gas,  and  hence  the  driving  field,  is  taken  to 
be  some  prescribed  function  of  time,  and  the  dynamical  equations  of  motion 
are  solved  subject  to  this  boundary  condition.  This  j:.«.s sumption  limits  the 
validity  of  the  analysis  to  the  time  during  which  the  prescribed  current  can 
be  maintained  experimentally  —  usually  not  the  entire  period  of  operation. 

It  also  involves  a  considerable  loss  of  detail,  such  as,  for  example,  the 
division  of  energy  at  arbitrary  times  between  magnetic  field  and  mechanical 
motion.  If  a  detailed  description  of  such  physical  quantities  is  desirable  foi 
design  or  diagnostic  purposes,  the  reaction  must  obviously  be  taken  into 
account. 

The  final  approximation,  that  of  idealizing  the  mechanical  properties  of 
the  gas,  consists  of  replacing  the  full  fluid  dynamic  equations  by  a  single 
equation  based  on  some  simplified  mechanical  model.  One  such  model  in 
current  use,  the  snowplow  model,  will  be  discussed  below.  Such  idealizations 
have  the  advantage  of  allowing  a  rather  simple  calculation  of  the  gross  dy¬ 
namics  of  the  gas  enveloped  by  the  shock.  However,  an  analysis  based  on 
this  kind  of  simplification  is  clearly  not  capable  of  providing  the  delailed 
description  of  the  flow  which  is  contained  in  a  solution  of  the  full  equations. 

The  latter,  unfortunately,  cannot  usually  be  solved  analytically. 

Of  the  assumptions  discussed  above,  the  first  two  are  essentially  un¬ 
avoidable  if  any  analytic  results  are  to  be  obtained.  Of  the  last  two,  at  least 
one  must  be  made  for  the  analysis  to  be  tractable.  If  both  (c)  and  (d)  are 
assumed,  the  analysis  may  become  particularly  simple  while  still  yielding 
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results  which  are  quite  adequate  for  many  purposes.  However,  if  quanti¬ 
tative  information  is  needed  concerning  the  interrelation  between  gas  and 
circuit,  the  reaction  of  the  gas  must  be  included.  The  desired  results  may 
then  be  obtained  with  somewhat  more  effort.  Finally,  there  are  special 
cases  where  similarity  solutions  of  the  full  fluid  dynamic  equations  can  be 
obtained  describing  the  flow  in  complete  detail.  To  obtain  these  solutions, 
it  is  again  necessary  to  replace  the  external  circuit  by  a  boundary  condition. 
Moreover,  the  form  of  the  boundary  condition  is  no  longer  arbitrary;  it  is 
determined  by  the  form  the  solution  must  take. 

Of  course,  in  addition  to  the  approximations  described  above,  which  are 
of  a  more  or  less  general  nature,  there  are  also  those  which  are  . specific  to 
a  particular  device. 

The  points  discussed  above  will  now  be  illustrated  by  specific  examples. 

II.  SNOWPLOW  THEORY  WITHOUT  CIRCUIT  REACTION. 

As  mentioned  in  the  Introduction,  the  greatest  analytical  simplification 
results  when  assumptions  (c)  and  (d)  are  both  made.  A  particularly  useful 
mechanical  simplification  is  the  so-called  snowplow  model,  which  was  first 
applied  by  Rosenbluth^  to  the  ordinary  pinch  effect.  It  is  assumed,  in  this 
model,  that  all  the  mass  swept  up  by  the  shock  is  compressed  into  an  in¬ 
finitely  thin  layer  immediately  behind  the  shock,  so  that  the  contact  front 
and  shock  are  the  same  interface.  The  motion  of  the  interface  is  determined 
from  the  principle  that  the  time  rate  of  change  of  momentum  of  the  accumu¬ 
lated  mass  is  equal  to  the  force  on  the  interface. 

It  is  interesting  to  note  that  snowplow  theory  can  be  derived  from  the  full 
gas-dynamic  equations  as  the  limit  shock  strength  °o,  specific  heat  ratio 
— *•!.  This  limit  is  called  Newtonian  theory^.  The  model  should  therefore 
be  a  valid  approximation  for  strong  shocks  in  gases  in  which  ionization  and 
dissociation  are  taking  place,  since  the  additional  degrees  of  freedom  pro¬ 
vided  by  these  processes  result  in  a  specific  heat  ratio  which  approaches 
unity.  It  is  clear  that  under  these  conditions  the  compression  (y  +  1)/(y  ~  1') 
becomes  very  large,  as  required  by  the  theory. 

In  the  devices  to  be  analyzed  below,  the  flow  is  approximately  one¬ 
dimensional;  that  is,  the  location  of  the  shock  at  a  time  t  is  describable  by  a 
single  coordinate,  X(t).  If  the  mass  of  gas  which  has  been  swept  up  by  the 
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shock  is  denoted  by  M(t),  the  basic  equation  of  snowplow  theory  can  be 
written 

-d-  (M  ^  )  =  F  .  (8) 

dt  dt 

where  F  is  the  force  on  the  interface.  This  equation  will  now  be  used  to 

3  4 

analyze  the  two  devices  known  as  inverse  pinch  and  Scylla  . 

A)  Inverse  pinch 

A  diagram  of  the  apparatus  for  this  device  is  shown  in  figure  1.  A  current 
is  passed  through  the  gas  and  returns  along  the  central  conducting  rod.  This 
current  produces  an  azimuthal  magnetic  field  which  pushes  the  gas  away  from 
the  rod,  leaving  behind  a  cylindrical  vacuum  region.  In  such  a  device,  the 
gas  can  be  pre-ionized,  and  an  axial  magnetic  field  produced  by  an  external 
solenoid  can  be  trapped  in  the  plasma  region.  The  resulting  shock  wave  will 
then  be  a  transverse  magnetohydrodynamic  one.  This  more  general  case 
will  in  fact  be  considered. 

The  additional  geometrical  idealization  will  be  made  of  replacing  the 
central  rod  by  a  line.  This  is  done  first  because  an  exact  solution  of  the  re¬ 
sulting  equation  can  then  be  obtained,  but  mainly  so  that  a  comparison  can  be 
made  with  the  similarity  solution  of  the  same  problem.  In  practice,  the 
solution  thus  obtained  should  become  a  good  approximation  to  the  flow  at 
distances  from  the  rod  of  several  times  its  radius. 

If  the  plasma  is  initially  at  a  uniform  pressure  p^  and  density  p^, 
and  in  a  uniform  axial  magnetic  field  of  strength  B^,  the  accumulated  mass 
per  unit  length  of  the  sheath  is 

M  =  TT  ,  (9) 


where  X  is  the  distance  of  the  sheath  from  the  axis,  while  the  net  outward 
pressure  on  the  sheath  is  (mks  units) 


P  " 


(10) 


The  first  term  on  the  right  hand  side  of  equation  (10)  is  the  magnetic  pressure 
of  the  aximuthal  field  produced  by  the  pinch  current  I,  while  the  second 
term  is  the  magnetic  pressure  of  the  external  field. 
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SCHEMATIC  DIAGRAM  OF  APPARATUS 
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Figure  1.  Diagram  of  apparatus  for  inverse  pinch 
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It  is  now  necessary  ho  make  some  assumption  concerning  the  form  of 
l{t).  Although  under  typical  experimental  conditions  the  current  is  usually 
sinusoidal,  the  ringing  time  of  the  circuit  can  be  made  sufficiently  long  that 
the  current  rises  linearly  over  a  large  portion  of  the  time  of  interest.  If  one 
then  assumes  a  linear  pinch  current,  I  =  wt,  and  introduce  as  parameters 
the  sound  speed  a^  =  (y  Pq/ Alfve^n  speed  , 

equation  (8)  with  F  =  ZitXp  becomes 


where 


jL. 

dt 


(X 


2^1,  ,,4  .x(b 

dt  '  o  X  o  Y  o 


C  =  (p.  I 
o  o 


2  2,  .  2 

03/471  p 


1/; 


(11) 


(12) 


is  a  characteristic  quantity  with  the  dimensions  of  a  speed.  The  solution  of 
equation  (12)  which  passes  through  X(o)  =  0  is 


X  =  kt, 
k2=  ^ 

K  -  4 


+ 


z 

Y 


+ 


r 


2,2 

a  ) 

o 


+  8  c 


(13) 


Thus,  according  to  the  theory,  the  front  moves  with  the  constant  speed 

dX  =  k. 
dt 


In  the  limit  of  strong  shocks,  where  the  snowplow  theory  should  give 
reasonably  accurate  results,  a^  /c^  <  1;  and  ^ 1;  equation 

(13)  becomes 


5” 


k  *  2 


(b^^  +  2_ 
o  ^  o 


(14) 


From  equation  (14)  and  the  definition  of  c^  the  scaling  laws  for  the  device  are 

obtained;  the  shock  speed  scales  as  the  square  root  of  the  rate  of  current  rise 

and  inversely  as  the  fourth  root  of  the  density.  A  comparison  of  the  snowplow 

12 

result  (for  a^  =  ^o  ~  experiments  of  Liepmann  and  Vlases  is  shown 

in  figure  2.  The  agreement  is  good  indeed,  verifying  the  scaling  law  with 
respect  to  both  current  rise  and  initial  density. 

In  addition  to  providing  the  scaling  laws  for  the  device,  the  theory  also 
provides  a  fair  idea  of  the  shock  speeds  to  be  expected  under  typical  ex¬ 
perimental  conditions.  For  example,  if  the  working  gas  is  deuterium  at 


26 


TN-61-29 


SUMMARY  OF  PRESSURE  PROBE  DATA 


Figure  2.  Comparison  of  snowplow  theory  of  inverse  pinch 
with  experiments  of  Liepmann  and  Vlases 
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100-|i,  initial  pressure  and  if  the  rate  of  current  rise  is  250  ka/n-sec,  the 
predicted  shock  speed  is  about  9  cm/p-sec,  again  in  good  agreement  with 
experiment. 

Finally,  although  the  theory  provides  no  information  about  the  detail  of 
the  flow,  the  internal  energy  in  the  gas  may,  within  the  framework  of  the 
theory,  still  be  calculated.  This  is  the  energy  which  is  available  to  heat, 
dissociate,  and  ionize  the  gas.  The  total  energy  per  unit  length,  E(t),  de¬ 
livered  to  the  gas  in  the  time  t  is  just  the  work  done  on  it  by  the  magnetic 
piston,  i.  e.  , 


X{t) 


X(t) 


E(t)  = 


F  dX  = 


(M  dX  )  jy 
dt  '  dt  ' 


X 


=  - 
dt 


M  dX 

dt^ 


(15) 


With  the  aid  of  equation  (8)  the  last  integral  can  be  transformed  to 


/ 


M  =  E(t) 

dt^ 


M 


■I 


'dt  ' 


dM 


so  that  finally 


M 


E(.l  =  \  M(^) 


(16) 


(17) 


The  first  term  on  the  right-haqd  side  of  equation  (17)  is  the  kinetic  energy  of 

the  gas;  the  second  term  is  the  internal  energy.  Internal  energy  is  thus 

1  Z 

acquired  at  the  rate  of  ^  (jt  ^  mass.  (The  dissipative  mechanism 

responsible  for  the  production  of  internal  energy  is  the  shock  wave.  )  The 
internal  energy  at  any  time  depends  on  the  history  of  the  motion  to  that  time, 
whereas  the  kinetic  energy  depends  only  on  the  velocity  at  that  instant.  When, 
as  in  the  case  under  consideration,  the  velocity  is  constant,  the  in- 
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tegral  in  equation  (17)  is  trivially  evaluated,  and  the  internal  energy  is 
evidently  just  equal  to  the  kinetic  energy.  In  such  cases,  then,  there  is 
exact  equipartition  between  kinetic  and  internal  energy. 

B)  Scylla 

In  this  device,  shown  in  figure  3,  the  external  circuit  drives  a  circum¬ 
ferential  current  around  the  outside  of  the  cylindrical  discharge  tube.  When 
the  switch  is  closed,  the  rising  current  induces  an  electric  field  which  drives 
a  surface  current  around  the  gas,  opposite  in  direction  to  the  primary  current. 
The  resulting  magnetic  field  is  axial,  and  drives  the  gas  radially  inward  from 
the  wall.  As  in  the  case  of  the  inverse  pinch,  such  a  device  can  be  used  to 
generate  magnetohydrodynamic  shock  waves  if  the  gas  has  been  pre-ionized 
and  an  axial  field  established  in  it  prior  to  the  discharge.  This  will  again  be 
the  case  considered. 

If  the  azimuthal  current  per  unit  length  is  denoted  by  i,  the  driving  axial 
field  is  uniform  between  the  current  sheet  and  the  outer  wall,  and  is  given  by 

=  fil  .  (18) 

and  the  snowplow  equation  of  motion  becomes 

ft  [.  =  -2.  X  ^  -pj  .  (19) 

where  X(t)  is  the  radius  of  the  current  sheet  and  X  is  the  radius  of  the  tube. 

o 

The  other  quantities  are  defined  as  before,  and  the  equation  is  to  be  solved 
with  the  boundary  condition  X(o)  =  X^  . 

An  analytic  solution  of  equation  (19)  can  be  obtained  if  it  is  assumed  that 
the  current  i(t)  rises  instantaneously  from  zero  to  a  constant  value  ij^.  (This 
can  be  approximated  fairly  well  in  practice  if  the  external  inductance  is  made 
sufficiently  small.  )  The  desired  solution  of  equation  (19)  is  then 

X  =  X  -  k, t  (20) 

o  1 
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where 

b  2  2  ll/2 

1  -  o  -  2-  (21) 

2  Y  2 

Cl  Cj 

and 


is  the  characteristic  speed  for  the  device. 

The  result  of  the  theory,  then,  is  that  the  gas  moves  radially  inward  with 
the  constant  speed  =  kj^.  It  is  interesting  to  note  that  the  speed  is  in¬ 
dependent  of  the  dimension  of  the  apparatus.  Moreover,  the  scaling  laws 
are  once  again  obtained;  the  shock  speed  scales  linearly  as  the  current  per 
unit  length  and  Inversely  as  the  square  root  of  the  density.  Once  again 
typical  shock  speeds  can  be  estimated;  for  deuterium  at  an  initial  pressure 
of  100-p  and  currents  of  about  10”  amps/ni,  shock  speeds  of  about  25-30 
cm/ p- sec  are  to  be  expected.  Finally,  the  internal  energy  may  be  calculated 
just  as  in  the  case  of  the  inverse  pinch;  here  again  there  is  equipartition  of 
energy. 

An  analysis  of  the  operation  of  either  of  the  above  devices  for  other  forms 
of  the  driving  current  in  general  requires  a  numerical  integration  of  the 
snowplow  equation.  The  same  is  true  if  the  solution  of  the  inverse  pinch  is 
desired  near  the  central  rod. 

III.  SNOWPLOW  THEORY  WITH  CIRCUIT  REACTION. 

As  mentioned  in  the  Introduction,  it  is  desirable  for  many  purposes  to 
have  a  detailed  description  of  the  interchange  of  energy  between  source  and 
device,  the  division  of  energy  between  magnetic  field  and  mechanical  motion, 
etc.  When  this  is  the  case,  the  analysis  must  be  modified  to  include  the  re¬ 
action  of  the  gas  back  on  the  energy  source.  This  modification  makes  it 
necessary  to  replace  the  single  equation  of  motion  of  the  last  section  by  two 
coupled  equations,  one  equation  for  the  electric  circuit  and  the  other  for  the 
mechanical  system.  A  detailed  analysis  of  this  type  has  been  carried  out 
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5 

for  the  ordinary  pinch  by  Killeen  and  Lippman. 

As  circuit  elements,  the  devices  in  question  are  essentially  variable 
inductances  which,  under  normal  operating  conditions,  can  be  considered  to 
be  connected  simply  in  series  with  the  energy  source.  ^  The  inductance  of 
the  device  at  any  instant  depends  on  the  location  of  the  plasma  boundary,  and 
it  is  this  dependence  which  couples  the  electric  circuit  and  the  mechanical 
system.  If  the  inductance  of  the  device  is  denoted  by  L^(t),  and  the  external 
source  is  taken  to  be  a  capacitor  of  capacitance  C  in  series  with  an  (un¬ 
avoidable)  inductance  L^,  the  circuit  equation  becomes 

+  ^  =  0  (23) 

where  Q(t)  is  the  charge  on  the  capacitor.  If  the  initial  pressure  of  the  gas 
is  neglected,  the  rate  at  which  mechanical  work  is  done  on  the  plasma  is 


dt 


(L 

( 


L.) 

1 


dQ 

dt 


W 


so  that  the  snowplow  equation  of  motion  (8)  becoqnes 


_sl_  /M  six, 
dt  '  dt  ' 


1 

2 


dLi 

“dX 


(dQ.)‘ 

'dt  ' 


(24) 


(25) 


The  system  is  thus  described  by  the  two  coupled  equations,  (23)  and  (25). 
In  these  equations,  dissipation  in  the  resistance  of  the  walls  of  the  device  and 
in  resistance  in  the  external  circuit  has  been  neglected.  Such  dissipation 
could  easily  be  included  but  in  practice  is  usually  small.  The  main  dissi¬ 
pative  mechanism,  the  shock  wave,  is  included  however. 

These  equations  will  now  be  applied  to  the  analysis  of  a  specific  device. 

MAST  (Magnetic  Annular  Shock  Tube^) 

This  device  (figure  4)  is  also  called  a  linear  plasma  accelerator.  The 
gas  is  contained  between  two  coaxial  electrodes  which  are  connected  to  the 
external  circuit.  When  the  switch  is  closed,  the  gas  breaks  down  near  the 
end  of  the  tube  and  a  radial  current  flows  in  a  thin  disk  between  the  inner  and 
outer  cylinders.  The  discharge  current  gives  rise  to  an  azimuthal  magnetic 
field  which  pushes  the  current  disk  down  the  tube.  The  disk  sweeps  up  the 
gas  ahead  of  it,  leaving  vacuum  behind. 
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LINEAR  PLASMA  ACCELERATOR 


Li  =  ix  Mv/X 

»mX 


Figure  4.  Diagram  of  apparatus 
for  linear  plasma  accelerator 
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Other  devices  for  the  linear  acceleration  of  plasma,  in  which  the  mass  of 

accelerated  gas  could  be  considered  constant,  have  been  analyzed  by  various 
7  8  9 

authors  '  '  along  these  same  lines.  These  analyses  differ  from  that  pre¬ 
sented  below  mainly  in  that  they  do  not  contain  the  mechanism  of  shock  dissi¬ 
pation.  One  treatment,  however,  ^  does  take  into  account  electrical  resis- 

9 

tance.  The  device  considered  here  has  been  discussed  by  Dattner  whose 
treatment,  however,  is  somewhat  inconsistent  in  that  it  is  assumed  that  the 
internal  energy  of  the  gas  is  always  some  fixed  fraction  of  the  kinetic  energy. 
It  is  clear  from  equation  (17)  that  this  is  not  generally  the  case. 

If  X(t)  denotes  the  distance  of  the  current  disk  from  the  end  of  the  tube, 
the  inductance  of  the  device  is 


where 


L.(t)  =  jt  X(t) 


(26) 

(27) 


is  the  inductance  per  unit  length,  d^  and  d^,  being  the  inner  and  outer  dia¬ 
meters,  respectively.  Furthermore, 

M(t)  =  p  A  X{t)  =  m  X(t)  (28) 

o 


where  A  is  the  area  of  the  annulus,  so  that  m  as  defined  by  equation  (28)  is 
the  mass  per  unit  length  of  the  tube.  If  the  initial  charge  on  the  capacitor  is 
Q^,  the  equations  of  motion  must  now  be  solved  subject  to  the  boundary 
conditions. 


Q(0)  =  Qq.  Q(0)  =  O.  X(0)  =  O,  X(0)  =  O  (29) 

A  study  of  the  problem  in  dimensionless  variables  shows  that  the 
operation  of  the  device  can  be  described  in  terms  of  a  single  parameter. 
If  one  introduces  as  dimensionless  variables 


(length) 

(time) 


(30) 
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q  = 


Q  (charge) 


the  equations  of  motion  become 


-d- 

dT 


[(1  +  x)  +  q  = 


dx  dx 

and  the  initial  conditions  are 


q(OI  =  -tS 


Q. 


y^2m 


=  q. 


(31) 

(32) 

(33) 


<1  (0)  =  x(0)  =  it  (0)  =  0 

The  equations  and  boundary  conditions  thus  contain  the  single  parameter  q^, 
the  ratio  of  the  initial  charge  on  the  capacitor  to  the  characteristic  charge  for 
the  system.  The  time-dependence  of  all  nondimensional  variables  is  therefore 
determined,  once  the  value  of  this  single  quantity  has  been  specified. 

A  general  analytic  discussion  of  the  system  (31)  and  (32)  presents  formid¬ 
able  difficulties.  However,  some  general  remarks  can  be  made.  For  very 
small  times,  x  <<  1,  that  is,  the  inductance  of  the  device  is  small  compared  to 
the  external  inductance.  The  external  circuit  is  thus  essentially  shorted 
through  the  device  and  starts  to  ring  at  its  natural  frequency.  From  the  re¬ 
sulting  driving  force,  the  initial  motion  of  the  gas  is  readily  calculated.  For 
X  <<  1,  the  solution  of  the  circuit  equation  is 

q  s  q^  cos  X ,  (34) 

and  with  this  as  the  driving  term,  the  solution  of  the  mechanical  equation  is 


/  2  2^/2 
X  =  ( X  -  sin  X ) 


(35) 


This  initial  phase  of  the  system  can  persist  for  a  considerable  time  if 
q^  is  made  sufficiently  small.  Since  a  small  q^  implies  a  weak  driving  force, 
the  condition  x  <<  1  can  then  be  satisfied  to  values  of  >>  1,  in  which  case 
the  velocity  of  the  front  ^  — *•  ^  .  If  x  becomes  too  large,  however,  this 

initial  description  breaks  down  and  a  period  of  strong  coupling  occurs. 
Ultimately,  however,  the  driving  force  (^  must  become  small;  the  gas 
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will  then  decelerate  while  conserving  its  momentum.  Thus,  in  the  final 
phase 

X  ~  kyr  (36) 

The  varying  inductance  now  "drives"  the  external  circuit,  the  equation  for 
which  becomes 


A. 

d  T 


- 

1 


dSL 

dx 


+  q  =  .0 


(37) 


The  asymptotic  solution  of  equation  (37)  for  large 


VM  A  0  j  0 


q 


X  sin 


(38) 


A  numerical  solution  of  equations  (31)  and  (3  \  has  been  carried  out  on 
an  IBM  7090  for  a  number  of  different  values  of  the  parameter  q^.  The 
results  are  shoWn  in  figures  5-8.  The  dashed  curve  is  the  asymptotic  solu¬ 
tion  for  X  <<  1  given  by  equations  (34)  and  (35).  Ii  is  evident  that,  for  the 
small  values  of  q^,  the  asymptotic  solution  is  a  yery  good  approximation 
indeed  to  the  actual  solution  over  the  entire  range  shown.  The  numerical 
solution  has  also  been  carried  to  values  of  x  >>  1;  the  results  verify  the 
validity  of  equations  (36)  and  (38)  in  this  limit. 

An  energy  integral  of  the  system  (31)  and  (32)  is  readily  obtained,  namely 


i 

4 


1 

2 


(38a) 

The  integral  in  equation  (38a)  is  the  internal  energy  in  nondimensional  form. 
The  internal  energy  can  be  calculated  from  this  equation  once  the  equations 
of  motion  have  been  integrated.  The  result  for  q^  =  1  is  shown  in  figure  9, 
where  for  comparison  the  kinetic  energy  has  also  been  plotted.  It  is  seen 
that  there  is  no  simple  relationship  in  this  case  between  internal  and  kinetic 
energy. 

The  predictions  of  the  theory  may  be  compared  with  the  experiments  of 

9 

Dattner.  These  experiments,  with  deuterium  as  the  working  gas,  covered 
a  range  of  pressures  of  200-1000  p  and  a  range  of  voltages  of  3-7  kv.  With 
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POSITION  OF  PLASMA  FRONT  AS  A  FUNCTION  OF  TIME 
FOR  VARIOUS  VALUES  OF  THE  INITIAL  CHARGE 
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VELOCITY  OF  PLASMA  FRONT  AS  A  FUNCTION  OF  TIME 
FOR  VARIOUS  VALUES  OF  THE  INITIAL  CHARGE 


T 


Figure  6.  Velocity  of  front  as  a  function  of 
time  for  various  values  of  q 

o 
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CHARGE  ON  CONDENSER  AS  A  FUNCTION  OF  TIME 
FOR  VARIOUS  VALUES  OF  THE  INITIAL  CHARGE 


Figure  7.  Charge  as  a  function  of  time 
for  various  values  of 
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CURRENT  AS  A  FUNCTION  OF  TIME  FOR  VARIOUS 
VALUES  OF  THE  INITIAL  CHARGE 


Figure  8.  Current  as  a  function  of  time  for  varioui 
for  various  values  of 

o 


40 


TN-61-29 


T 

Figure  9.  Energy  distribution  in  the  gas 
as  a  function  of  time  for  q  =  1 
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-6  -9 

circuit  parameters  of  C  =  12.  5  x  10  f  and  =  57,  4  x  10  h.  this  range 
of  pressures  and  voltages  corresponds  to  a  range  of  values  of  q  of  0.  8  -  4.  2. 
It  was  observed  that  the  velocity  of  the  gas  became  roughly  constant  after 
about  1/4  of  a  cycle,  with  the  velocities  attained  by  the  gas  ranging  from  6 
to  11  cm/fisec.  The  maximum  calculated  velocities  occur  at  about  the  proper 
time  and  range  from  4  to  14  cm/p  sec.  Thus,  there  Is  rough  agreement  be¬ 
tween  theory  and  experiment. 

One  of  the  sources  of  discrepancy  between  theory  and  experiment  lies  In 
the  assumption  that  the  flow  Is  one-dlmenslonal,  that  Is,  that  the  current 
Interface  Is  a  plane  disk.  The  magnetic  pressure  across  the  annulus  varies 
as  1/r^,  where  r  Is  the  radius,  and  therefore  decreases  from  the  Inner 
electrode  to  the  outer.  This  pressure  distribution  Is  Incompatible  with  the 
uniform  dynamic  pressure  which  exists  across  a  plane  disk.  It  Is  actually 
not  difficult  to  modify  the  theory  to  take  this  effect  Into  account.  The  shape 
of  the  Interface  can  be  calculated^^  by  requiring  that  the  magnetic  pressure 
across  the  Interface  Is  everywhere  balanced  by  the  dynamic  pressure  given 
by  snowplow  theory.  The  shape  calculated  In  this  manner  turns  out  to  be 
a  paraboloid  of  revolution.  In  agreement  with  experiment. 

IV.  SIMILARITY  SOLUTION  FOR  INVERSE  PINCH. 

As  mentioned  In  the  Introduction,  special  cases  exist  for  which  a 
similarity  solution  of  the  full  magnetohydrodynamlc  equations  may  be  ob¬ 
tained  describing  the  flow  In  complete  detail.  Of  course,  the  external 
circuit  must  be  replaced  by  an  appropriate  boundary  condition.  One  such 
case  Is  that  of  the  Inverse  pinch  with  just  those  Initial  and  boundary  conditions 
for  which  the  snowplow  solution  was  obtained  above,  namely  an  Initial  uni¬ 
form  axial  field  In  the  plasma  and  a  discharge  current  which  Increases 
linearly  with  time.  Details  and  results  of  the  calculation  appear  elsewhere. 

It  Is  sufficient  to  say  that  the  solution  Is  of  the  form 


p(x,  t)  =  p^o  (<?>) 
u(x,  t)  =  C*  U 
p  (x.  t)  =  p  ^C*^  P(<I>) 

B^Cx.  t)=p^c*^  M®) 


(39) 
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SPEED  OF  CONTACT  FRONT  AS  A  FUNCTION  OF 
INITIAL  PRESSURE  IN  GAS  FOR  DIFFERENT 
VALUES  OF  THE  RATIO  OF  INITIAL  MAGNETIC 
PRESSURE  TO  GAS  PRESSURE 


Figure  10.  Comparison  of  speed  predicted  by  snowplow  theory  with 
speed  of  contact  front  obtained  from  a  similarity  solution  of  the  full  equations 
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where  o ,  U,  P,  and  p  are  nondimen sional  functions  of  a  single  nondimension- 
al  variable  <1'.  Two  exact  integrals  of  the  four  magnetohydrodynamic  equa¬ 
tions  are  obtained,  and  the  remaining  two  equations  integrated  numerically. 
The  full  detail  of  the  flow  between  the  shock  and  the  contact  front  is  thus  ob¬ 
tained. 

Of  particular  interest  is  the  result  shown  in  figure  10.  Here  the  speed 
of  the  contact  front  obtained  from  the  solution  of  the  full  equations  is  com¬ 
pared  with  tbe  speed  predicted  by  the  snowplow  theory,  equation  (13).  It  is 
seen  that  there  is  remarkably  good  agreement  between  the  snowplow  solution 
and  the  similarity  solution  over  the  entire  range  of  shock  strength.  The 
reason  for  this  uniformly  good  agreement  is  that,  independent  of  shock 
strength  and  specific  heat  ratio,  most  of  the  momentum  in  the  flow  is  carried 
by  the  fluid  very  close  to  the  contact  front  where  the  density  and  velocity  are 
greatest.  Snowplow  theory,  being  simply  a  statement  of  momentum  balance, 
predicts  rather  accurately  the  motion  of  this  portion  of  the  fluid.  This 
accounts  for  the  success  of  the  theory  noted  in  the  preceding  sections. 
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NUMERICAL  APPROXIMATIONS  FOR  WEAK  SOLUTION  OF  MIXED 
INITIAL-BOUNDARY  VALUE  PROBLEMS  OF  FLUID  FLOWS 
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Mr.  M.  Halem 

Republic  Aviation  Corporation 


1.  INTRODUCTION. 

This  paper  will  present  a  brief  survey  of  some  of  the  problems  that  are 
under  Investigation  by  the  Plasma  Propulsion  Lab,  in  conjunction  with  the 
Digital  Computing  Division  at  Republic  Aviation.  Each  of  the  problems 
represent  a  special  case  of  a  more  general  system  of  conservative  hyper¬ 
bolic  equations  with  moving  boundaries. 

In  support  of  the  work  being  carried  out  in  developing  a  pulsed  plasma 
accelerator,  a  number  of  studies  have  been  made  of  the  pinch  process  using 
various  models  with  varying  degrees  of  sophistication.  In  the  present  paper 
emphasis  will  be  placed  on  the  analytical  solution  obtained  by  Chu^  and  th6 
numerical  gasdynamic  treatment  using  the  technique  of  the  pseudo  viscosity 
term  of  Lax^  lor  application  to  an  IBM  7090.  The  analytical  approach  is 
based  upon  the  assumption  that  one  shock,  and  only  one,  exists.  The  numer¬ 
ical  approach  makes  no  such  restriction.  It  is  conceivable  in  the  general 
case  that  more  than  one  could  exist,  resulting  in  shock  collisions  and  violation 
of  the  constant  entropy  condition  for  piston  problems.  On  using  the  numer¬ 
ical  technique,  however,  if  more  than  one  shock  does  arise  they  are  im¬ 
mediately  recognizable  as  areas  across  which  rapid  changes  in  the  variables 
take  place,  and  are  automatically  accounted  for  in  the  analysis.  This  tech¬ 
nique  is  therefore  a  very  powerful  one,  especially  in  cases  where  the  actual 
physical  processes  are  not  completely  understood.  This  same  method  has 
been  applied  to  the  study  of  hydromagnetic  shocks  generated  by  the  relative 
motions  of  fluids  in  the  presence  of  magnetic  fields.  Here  again  the  shocks 
are  readily  identified  and  can  be  categorized  into  the  fast,  slow,  or  Alfven 
variety,  without  the  necessity  of  prior  shock  fitting.  The  technique  is 
applicable  to  fluid  flows  for  arbitrary  dimensions. 

The  analytical  approach  is  presently  being  extended  to  the  compression 
of  a  sphere  of  plasma  by  an  external  magnetic  field,  where  the  field  intensity 
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is  a  function  of  the  plasma  shape  at  any  instant.  This  problem  is  being 
studied  in  support  of  project  Cantaloupe  under  contract  with  AFSWC. 

A  further  problem  that  fits  into  the  general  numerical  procedure  is  the 
Stefan  problem  of  heat  conduction,  which  is  extremely  important  in  reentry 
studies.  The  program  that  has  been  developed  makes  it  possible  to  treat 
arbitrary  geometries  in  one,  two,  or  three  dimensions,  with  layers  of 
different  materials  and  moving  boundaries.  A  similar  problem  also  arises 
in  the  studies  of  the  initial  ionization  processes  and  skin  formations  in  the 

3 

pinch,  as  was  described  by  Killeen,  Gibson  and  Colgate  . 

2.  1  HYPERBOLIC  SYSTEMS  OF  CONSERVATION  FORM. 

Consider  the  hyperbolic  systems  of  equations  from  gasdynamics  and 
magnetogasdynamics  which  have  the  following  form  governing  the  basic  flow 
equations: 

Ut-V-F^-B  =  0 

where  U  is  a  column  vector  of  the  unknown  functions,  F  and  B  are 
vector  valued  functions  of  x^,  t  and  U  .  where  x^  represents  the  spatial  co¬ 
ordinates  X,  y,  z  and  the  symbol  V  is  the  customary  del  operator.  The 
initial  condition  U  (x,  0)  ■  <{)(x),  is  assumed  to  satisfy  equation  (2.  1).  Such 
systems  of  equations  are  said  to  be  of  conservative  form. 

To  illustrate  such  systems,  consider  the  nonsteady,  non-isentropic, 
dimensionless  ideal  gas  flow  equations  in  Eulerian  coordinates  which  are 
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where 


and  where  ^  O  ^ 

u ,  r  y  are  dimensional  quantities 
and  <3  m  0,  1,2  for  the  one-dimensional,  axisymmetric,  spherically  sym¬ 
metric  cases  respectively. 

The  familiar  Lundquist  equations  (for  infinitely  conducting,  nondissipative, 
inviscid  ideal  gas)  in  dimensionless  form  are; 


By 

di 


a 

ai 


(2.  1.2) 
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=0 


-a-Tio!: 

3t 


^(mx+B^rfiY 

«) 


where 


ia^lf 

Ko-i 


3 


\i  =.]  ^ 

.U^Y. 


and  U  is  the  initial  value 
yo 

of  U  .  • 

y 


The  equations  for  a  steady,  two-dimensional,  nonviscous,  adiabatic, 
rotational  flow  are:  . 

V  (^uJ)'  o 


^(a,^v)+^(vpv)_^p  =0 

^  ^  pi  O 
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where  p,  p ,  u  are  the  unknowns. 

These  particular  systems  are  only  a  sample  of  those  which  lend  them- 
selves  to  a  general  method  of  numerical  solution  which  is  directly  applicable 
to  the  actual  form  of  the  physical  problems.  The  occurrence  of  these  prob¬ 
lems  in  nature  usually  involves  at  least  one  moving  boundary  and  more  often 
two  free  boundaries.  Almost  any  attempt  at  their  solutions  has  led  to  simpli¬ 
fying  assumptions  which  usually  are  not  physically  realizable. 

For  such  systems  a  complete  mathematical  theory  exists,  since  they  are 
all  of  hyperbolic  type.  However,  the  existence  and  uniqueness  of  solutions 
of  such  systems  are  all  of  a  local  nature.  The  physical  nature  of  the  problems, 
because  of  neglecting  such  notions  as  viscosity  and  heat  conduction,  gives  rise 
to  the  mathematical  notion  of  "shocks"  or  discontinuities.  These  discontin¬ 
uities  in  the  solutions  do,  however,  have  certain  algebraic  properties  derived 
from  the  conservation  principles.  It  is  these  curves  which  become  the  free 
boundaries  in  many  important  physical  problems  and  they  have  led  to  a  host 
of  analytical  and  numerical  techniques  for  effecting  their  solutions. 

2.  2  THE  "LAX  DIFFERENCE  SCHEME". 

The  method  adopted  for  the  numerical  solution  of  such  a  class  of  free 
boundary  problems  is  essentially  due  to  Lax^.  The  notion  of  weak  solution 
is  introduced,  since  one  is  interested  in  solutions  in  the  large  which  need  not 
have  continuous  first  derivatives.  A  solution  is  defined  as  weak  if  it  satisfies 

|j Kv  +FV'W-Vi6]  jv/(x,d)0U)cix  =  0  (2.21) 

where  W  is  a  test  function  with  continuous  first  derivatives  and  is  zero  out¬ 
side  some  interval.  It  can  be  shown  that  for  two  genuine  solutions  on  two 
adjacent  domains  separated  by  a  smooth  curve,  where  the  slope  of  the  curve 
satisfies  the  Rankine-Hugcniot  shock  conditions  with  respect  tc  these  two 
solutions,  the  smooth  curve  forms  a  weak  solution  over  the  entire  domain. 

As  so  often  happens  when  one  generalizes  the  notion  of  solutions  in  mathe¬ 
matics,  a  penalty  is  imposed.  In  this  case  the  initial  values  for  conservative 
systems  do  not  determine  a  unique  weak  solution.  However,  an  additional 
principle  asserted  by  Lax  to  determine  the  unique  solution  is  that  "weak 
solutions  for  fluid  problems  are  limits  of  viscous  flows". 
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Using  this  principle  one  can  enlarge  equations  (2.  21)  to  a  nonlinear 
parabolic  system 

Ut  +  VF +B  =  (2 


For  a  fixed  initial  value  of  U(x,  o),  the  solution  exists  for  a  range  of  t 

independent  of  X.  .  If  corresponding  solutions  (x,  t)  converge  boundedly 

to  a  limit  U(x,  t)  as  X  -•  0.  then  this  limit  will  be  a  weak  solution. 

The  difference  scheme  is  derived  by  choosing  .  ( A r)^  and  using  a 

^  "  2At 

forward  difference  approximation  with  respect  to  t  and  centered  difference 
approximations  with  respect  to  x. 

This  difference  method  gives  rise  to  the  following  recipe  for  differencing 
the  equations. 


Replace  d  U  by 
t 


At 


■a(x/+A4:')  -  U  (x  +  U 


and 


d  F 


d  X 


by 


J 


Ax 


f6j  -  F  Dj(x-AX.t^l 


(2.2.3) 


(2.  2.  4) 


This  difference  scheme  when  viewed  as  an  approximation  to  equation  (2.  2.  2) 
satisfies  the  stability  condition  for  parabolic  equations; 


i.  e. 


z 


The  difference  equations  can  also  be  shown  to  be  stable  for  the  hyperbolic 
systems. 

Since  it  is  desired  to  select  the  unique  solution  in  nature,  it  is  con¬ 
jectured  that  one  should  try  to  satisfy  the  stability  criteria  of  Courant- 

Friedricks-Lewy  for  hyperbolic  systems  (i.e.  ,  <  _ L _  ). 

Ax  u  +  c 

It  is  our  aim  to  extend  this  method  to  mixed  initial  value  problems.  For 
free  and  fixed  boundary  value  problems  one  usually  has  integral  forms  of  the 
conservative  systems  determining  the  behavior  of  the' boundaries.  The 
system  of  integral  equations  take  the  following  form 
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f 

Vft) 


(Z.i.S) 


These  integral  forms  lead  to  equations  (2.1.  l)by  application  of  Gauss' 
Theorem.  However,  the  equations  being  solved  are  of  the  form  (2.2.2).  Thus, 
equations  (2.2.  5) should  be  modified  if  they  are  to  lead  to  equation  (2.  2.  2).  A 
natural  generalization  of  the  integral  forms  leading  to  this  system  takes  the 
form 


(  UtcfY*  ( 


FaciS  = 

S(t) 


AV\J^.(5S 

'5(0 


(2.  2.  6) 


From  this  form  it  is  possible  to  obtain  equations  governing  the  motion  of 
the  free  boundaries. 

(^estions  as  to  the  existence  of  solutions  of  equations  (2.  2.  2)subject  to 
conditions  (2.  2.  5)  have  been  only  partially  answered  for  specific  problems. 

That  such  solutions  of  the  equations  (2.2.  2)  with  boundary  conditions 
(2.  2.  ^converge  to  solutions  of  equations  (2.  1. 1)  subject  to  boundary  conditions 
(2.2.5)can  only  be  conjectured  at  this  moment.  However,  for  many  of  the 
cases  tried  so  far,  the  method  seems  to  give  promising  results. 

3,  APPUCATIONS  TO  PHYSICAL  PROBLEMS. 

The  general  difference  scheme  just  described  will  be  applied  to  various 
problems  which  have  already  been  solved  by  special  techniques.  A  summary 
of  tho  technique  by  which  analytic  solutions  have  been  obtained  will  be  de¬ 
scribed,  and  then  the  application  of  the  generalized  Lax  method  to  the  same 
problem  will  be  presented.  The  following  selection  of  problems  includes  some 
of  those  that  have  been  of  theoretical  concern  to  the  Plasma  Propulsion  Lab 
and  of  numerical  interest  to  the  Digital  Computing  Division  of  Republic 
Aviation  Corporation. 

3. 1  STEFAN' S  PROB  LEM. 

The  first  type  of  application  to  be  considered  includes  those  physical 
problems  governed  by  parabolic  equations  with  moving  boundary  conditions. 
Missile  reentry  problems  for  bodies  with  ablating  materials,  melting  of  ice, 
recrystallization  of  metals,  are  some  examples  of  problems  now  considered 
as  Stefan's  problem.  The  problem  also  arises  in  studies  of  the  initial 

3 

ioniaation  processes  in  a  plasma  .  Consider  the  simplest  form  of  the  Stefan 
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problem  uescribed  by  the  following  equations; 


A 


ax2- 


(3.1.1) 


for  O  <  X  <1  t  (0  where  b  (t)  is  the  position  of  the  boundary  subject  to 
the  following  initial  and  boundary  conditions 


t  >0 


vfx.o)  =  X 


(3.  1.2) 


V  =  1 


lx  ^ 

b(o)  =  1 


4 

This  is  almost  the  form  considered  by  Douglas  and  Gallie  for  which  a 
numerical  difference  scheme  is  shown  to  converge  to  a  solution  whose  ex- 

5 

istence  and  uniqueness  have  already  been  established  by  Evans  Theirs 
is  a  series  method  that  resolves  the  solution.  However,  the  coefficients 
of  this  series  are  difficult  to  obtain  and  no  radius  of  convergence  is  available. 

Applying  the  recipe  for  the  equation  of  form  (2,  2.  2)to  the  equations  (3.  1.  1) 
where  U  ■  v  and  F  -  0  leads  to  the  following  difference  equation  at  some 
point  X,  t  in  the  interior  of  the  domain  Q  —  ^  ^ 

v(x,-t+4t)  =  2vrx,t)  +  v(K-AX,i')7  ,3  J  3, 


Additional  equations  are  now  needed  for  the  boundary  points  [b  (t).  t]  . 
Equations  (2.  2,  6)applied  to  (3.  1.  l)lead  to  the  following  equations 

f  'bit) 


dt 


+  A  ■= 


ix 

at 


cSx 


(3.  1.4) 
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Since  v  Lb(t),  tj  ii  known  for  all  t,  (3. 1. 4)givee  an  equation  lor  b(t). 
Approximating  the  integral  by  a  trapeaoidal  rule  and  replacing  time  de¬ 
rivatives  of  the  quantities  by  forward  differences  lead  to  the  following 


■bft  1- At)  =  b(t)  +  At  fi  - X 

whtre  N  AX  =  bW 


(3.1.5) 


The  general  recipe  for  the  solution  of  free  boundary  value  problems  for 
hyperbolic  systems  obtained  from  the  integral  conservation  laws  was  an 
outgrowth  of  this  problem.  The  effectiveness  of  this  approach  for  a  single 
nonlinear  parabolic  equation,  with  boundary  conditions  general  enough  to  in¬ 
clude  radiation  terms  and  ablation,  is  described  in  an  RAC  report  by  Pines, 
Halem  and  Broder.  ^ 


It  was  natural  to  try  to  extend  this  approach  to  nonlinear  parabolic  sys¬ 
tems.  Other  numerical  schemes  for  such  systems  are  presently  being  in¬ 
vestigated,  and  as  yet  many  open  questions  remain  to  be  answered. 


3.  2  HYDROMAGNETIC  SHOCKS 

The  next  class  of  problems  that  the  general  numerical  procedure  will  be 
applied  to  are  the  problems  leading  to  magnetic  shocks  caused  by  relative 
motions.  Consider  a  conducting  pole  face  of  an  electromagnet  adjacent  to  an 
Infinitely  conducting  fluid  with  a  magnetic  field  normal  to  the  face  and  where 
no  relative  motion  between  magnet  and  fluid  is  assumed.  (Figure  1. ) 

Suppose  the  magnet  is  suddenly  set  in  motion  with  a  velocity  '  ^ 

motion  will  result  because  of  the  presence  of  the  x-component  of  the  magnetic 
field.  The  magnetic  lines  of  force  can  be  regarded  as  fixed  between  magnet 
and  fluid.  The  motion  of  the  magnet  pulls  down  the  lines  of  magnetic  inten¬ 
sity  and  the  fluid  in  a  thin  layer  adjacent  to  the  pole.  This  produces  a  com¬ 
ponent  of  the  magnetic  field  in  the  layer.  The  variation  of  B^  with  x  in 
le  layer  results  in  a  current  density  in  the  z  direction  of  magnitude 


dH 


Thus  the  fluid  experiences  a  magnetic  force  J  xlS  which  is  in 


the  X  direction  of  magnitude  B 


d  H 


This  force  produces  a  motion 


in  the  layer  which  then  propagates  out.  One  can  consider  the  motion  of  the 
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Figure  1. 


|\/o  OF  CVCLES 


± 
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Figure  Z. 
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electromagnet  in  the  x  direction,  which  leads  to  a  pure  gas  dynamic  problem, 
or  in  some  arbitrary  direction  leading  to  more  complicated  shock  and  wave 
patterns. 

The  resultant  flows  will  give  rise  to  fast,  slow,  and  Alfven  hydromagnetic 
shocks  as  well  as  various  wave  motions.  The  particular  problem  for  motion 
in  the  Y  direction  was  firstTormulated  by  K,  O.  Friedrichs  and  leads  to 
resolution  of  shear  flow  discontinuities.  J.  Bazer  was  the  first  to  obtain 

7 

numerical  solutions  of  these  problems  .  The  method  used  was  that  of  piecing 
together  fast  shock  and  simple  wave  solutions. 


The  general  approach  described  in  Section  2,  when  applied  to  the  Lund- 
quist  equations  (2.  1.2)  leads  to  the  following  system  of  difference  equations. 
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(3.  2.  Dcon't. 


Select  a  coordinate  lyetem  that  moves  with  the  mean  velocity  of  the 
thin  layer.  The  initial  vector  is  U  ^  p„0  ’^‘^0 

find  U  C,  >'^x X>0.  The  numerical  solution  with 
this  difference  method  gave  very  good  agreement  with  the  numerical  an¬ 
alytical  solution  of  Baser  (figure  2).  For  better  agreement  in  the  neighbor- 
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hood  of  the  shock,  larger  values  of  —  are  needed  whereas  for  the  steady 
state  region,  smaller  mesh  spacings  should  be  used. 

This  approach  is  now  being  used  in  resolving  other  hydromagnetic  shock 
flows  and  seems  very  promising. 

The  method  gives  numerical  answers  which  are  sufficiently  accurate  to 
be  used  directly  in  shock  fitting  techniques.  Quite  often  the  knowledge  of  the 
strength  of  some  parameters  across  magnetic  shocks  of  arbitrary  types,  as 

g 

well  as  the  relations  which  are  known  ,  enables  one  to  actually  compose  the 
flows  from  the  pieces  directly.  These  techniques  are  now  being  developed 
i>r  arbitrary  motions  of  the  magnetic  poles  moving  into  infinitely  conducting 
plasmas. 


3.  3  MAGNETIC  PISTON  PROBLEM  (PINCH  EFFECT) 


Consider  a  pair  of  circular  electrodes  in  an  electrical  circuit  containing 
capacitive,  resistive,  and  inductive  elements.  When  a  sufficiently  high 
voltage  is  applied  across  electrodes,  the  gas  between  them  ionizes  and  a 
current  sheath  confined  to  a  thin  cylindrical  shell  at  the  edge  of  the  electrodes 
is  produced  enclosing  column  of  gas  between  electrodes.  The  magnetic  field 
produced  by  current  sheath  then  interacts  with  the  current  to  produce  a 
adially  inward  Lorentz  force,  J  x  B,  and  compresses  the  enclosed  gas. 

If  the  ionized  gas  has  high  conductivity,  a  current  sheath  of  approximately 
zero  thickness  is  formed,  which  enables  one  to  think  of  this  sheath  as  a 
magnetic  piston.  Since  the  magnetic  pressure  and  gas  pressure  should  be 
continuous  across  the  infinitely  thin  current  sheath,  we  obtain  gas  pressure 


—  "in-  ^  ^  where  Pg  i»  static  gas  pressure  behind  the  piston 


'9  8  7‘r 


which  varies  in  time.  It  will  be  assumed  initially  that  I  (t)  is  given  from 
external  circuit  characteristics  and  the  time  varying  inductance  due  to  move¬ 
ment  of  sheath  is  neglected.  Assuming  no  leakage  through  the  magnetic 
piston  (or  current  sheath)  equations  (2. 1.1)  govern  the  flow.  Prescribing  p^ 
on  the  curve  OP  in  (figure  3)  and  the  initial  state  on  OA  is  a  well  poted 
problem. 

BS  represents  the  shock  curve  arising  from  the  compression  by  the 
piston.  The  positions  of  the  curves  OP  and  BS  as  well  as  the  flow  between 
them  is  to  be  determined. 


58 


TN-61-29 


Several  approximate  methods  for  solving  this  problem  have  been  de¬ 
veloped,  but  a  complete  solution  to  the  physical  problem  is  still  open,  A 
widely  used  technique  for  approximating  the  solution  is  the  snowplow  model. 
This  thoery  assumes  that  all  the  mass  swept  up  by  the  shock  is  compressed 
Into  a  thin  layer  so  that  the  piston  and  shock  curves  can  be  treated  coinci¬ 
dentally.  The  motion  of  this^layer  is  governed  by  the  momentum  equation. 
This  model  gives  valid  approximations  when  the  shock  is  sufficiently  strong 
and  the  specific  heat  ratio  of  the  gas  is  close  to  unity. 

Another  approach  is  that  described  in  the  paper  by  Drs.  Greifinger  and 

9 

Cole  which  makes  use  of  similarity  solutions.  Assuming  the  solution  has 
conical  similarity,  and  introducing  a  coordinate  system  employing  the  stream 
function  and  a  certain  mass  ratio  parameter,  one  obtains  two  coupled  ordinary 
differential  equations. 

3.  3A  An  Analytical  Formulation  of  the  Pinch  Problem. 

Employing  a  somewhat  different  approach  by  B.  T.  Chu^  it  is  also 
possible  to  reduce  the  partial  differential  equations  governing  the  flow  to  a 
coupled  system  of  ordinary  differential  equations.  By  labeling  the  particles 
as  when  they  are  traversed  by  the  shock  (figure  4) 
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Figure  4. 

aLagrangian  formulation  of  the  equations  for  cylindrical  geometry  is 
obtained  in  the  form 


/• 

I 


where  the  subscript  S  indicates  that  the  quantity  to  which  it  is  affixed  is 
evaluated  just  behind  the  shock.  For  a  strong  shock  the  density  ratio 

The  boundary  condition  to  be  satisfied,  obtained  from  the  usual  shock  con¬ 
ditions  and  magnetic  piston  in  dimensionless  form,  are 
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r^p  -  L*(t) 

Q  f  =0 

p  “  Lbf 

®  Y  =t 

I 

p  ~  ^  i 

®  v|)  ^-b 

(3.  3.  2) 

*=  2  n 

(g)  ^ 

r- 1 

@  ^  *  t  -  0 

where  U  is  the  shock  speed  (Note:  ~ 

Assuming  strong  shocks  one  can  write  the  solution  in  a  power  series  e 

p  =  p‘°^  +  ep^'^ 

+  G  p  • 

r  =  +  er"' 

(3.3.3) 

t  .  ..  ) 

u  =  11“’’+  euf'^  -^e"ui*’+  ■  ■ 

u.„-  +  £U,“  +eVf  + .  . . 

Substituting  into  the  equations  (3.3.  2)and  collecting  like  terms  of  e,  the 
lowest  order  terms  are  found  to  satisfy  the  differential  equation 


(3.  3.4) 


WHEUe  r‘°^  =  fU)  ANO 
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This  equation  is  identical  to  that  obtained  from  the  snowplow  model. 

Since  the  initial  condition  of  f»l  at  t«0  is  singular,  the  solutions  of  these 
equations  are  determined  by  the  single  condition  that  f"  (0)  <0.  When  I(t)  ■ 
sin  cC  t  one  can  derive  a  power  series  approximation  for  f  (t)  which  for 
small  oCt  is 

Kt)  ’  I  "  ^  (3-3-5) 

from  which  one  obtains  f"  (0)  as  a  starting  condition. 

Making  the  substitutions  and  equating  first  order  terms  then  leads  to  the 
following  integrodifferential  equation 

*^0  -Jo 

F(ni)--rw7[r^.r(fu)-j{va 

GCW  =  -  f  1  K'f'.O.iH'l 

a*'-'*  I'-F  +  rFCf,t)c54> 
al"-  flv) 

At  t  s  0,  f  >  1  is  a  singular  point,  and  the  solution  is  determined  from  the 
condition  f^  (o)  m  o. 
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To  obtain  initial  values  one  again  obtains  a  power  series  representation 
of  (t)  for  small  0C  t  as 

The  numerical  technique  used  to  solve  the  coupled  system  of  differential 
equations  was  a  general  self- starting,  Cowell  (double  sum)  second  order  in¬ 
tegration  scheme  with  eighth  order  accuracy.  The  routine  has  the  option  of 
using  a  fixed  step  or  variable  step  size  for  better  corrections.  This  method 
is  more  rapid  and  has  greater  accuracy  than  Runge-Kutta  schemes.  The 
integrals  are  evaluated  by  Simpsen's  rule  and  the  partial  differential  operator 
on  the  integrals  by  a  backward  second  order  difference  scheme.  The  power 
series  solution  for  small  oCt  is  used  to  allow  for  a  sufficient  number  of  points 
to  build  up  good  integral  and  differential  approximations.  The  integration 
scheme  is  used  with  a  sufficiently  small  fixed  step  size  to  avoid  exceeding 
the  memory  storage  requirements  of  the  7090  which  might  occur  if  variable 
time  steps  were  used.  The  size  of  the  steps  was  determined  by  solving 
special  cases  of  these  equations  with  the  variable  time  step  option.  The  so¬ 
lutions  give  the  snowplow  results  as  the  zeroth  order  approximation,  while 
the  first  order  result  gives  the  separation  distance  betwen  the  piston  and  the 
shock  (figure  5.  ).  This  distance  is  relatively  small,  so  that  the  dynamics  are 
v«ll  satisfied  by  the  snowplow.  However,  the  first  order  solution  gives  the 
distribution  of  density,  pressure,  temperature  and  velocity  in  this  small 
piston- shock  space.  Because  of  the  assumption  of  strong  shocks,  the  density 
at  the  piston  is  infinite  initially,  and  remains  so  for  the  rest  of  the  motion. 

The  basic  assumption  of  such  analytical  solutions  is  that  only  one  shock 
is  formed  during  the  piston  acceleration.  It  is  conceivable  that  a  number  of 
shocks  could  develop  in  the  general  case,  and  under  these  conditions  re¬ 
flected  waves  could  change  the  constant  entropy  condition  at  the  piston.  It 
is  in  this  respect  that  the  Lax  procedure  represents  an  extremely  powerful 
tool,  since  high  gradients  throughout  the  flow  field  will  become  immediately 
evident. 
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Figure  5.  Pinch  dynamics  on  a  constant  area  channel 
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3.  4  THE  LAX  DIFFERENCE  APPROXIMATIONS  FOR  THE  PISTON 
PROBLEM. 

F or  the  present  purposes  the  difference  equations  will  be  specialized  to 
treat  one-dimensional  flow. 

Consider  a  subdivision  C  of  the  domain  D  of  figure  3  and  denote  the  mesh 
points  Y  *  and  "t  A^~f.  where  i  and  j  are  non-negative 

integers.  Let  U  (r,  t)  at  these  points  be  denoted  as  U..  Using  the  recipe  de¬ 
scribed  in  (2.2.3,  2,  4)  the  following  system  of  difference  equations  over  the 
interior  points  with  ^  andCT’^O  can  be  derived  from 

equations  (2.  1.  l)as 


(3.4.1) 


-  a- 

Z  p-.  j 
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The  time  centered  difference  approximation  over  unequally  spaced 
distance  intervals  Ar  and  ocLr  at  the  mesh  points  adjacent  to  the  piston 
path  (figure  6)  can  be  derived  from  the  Taylor  series  expansion  in  both  di¬ 
rections;  and  eliminating  the  second  order  terms,  one  obtains  for  a  function 


<P  (r,  t) 


and  similarly 


(3.4.2) 


oC(oC 


(3.  4.  3) 


oc(qc^  -I) 

cC  ^  i)(^  2.)  Lr 


(l->) 


(t-o 


Apply  these  approximations  to  (2.1.  Done  obtains  the  following  system  for 
the  points  adjacent  to  the  piston  path: 
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3.  5  NUMERICAL  APPROXIMATIONS  AT  THE  BOUNDARY. 

The  fundamental  disturbances  of  the  fluid  arise  and  propagate  from  the 
moving  boundary  and  present  the  inherent  difficulties  of  the  problem. 

Various  numerical  techniques  were  attempted  and  the  appropriate  approxima¬ 
tions  will  be  discussed. 


3.5A 


The  generalized  system  of  integral  equations  to  be  used  is: 


j-r(*)  (3.5A1) 

ca)  cir  =  I  ^  dr  \  (A  ?r 

*^0  -'o 
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The  additional  boundary  conditions  arise  from  the  equation  of  state  and 
the  equation  for  the  magnetic  pressure. 

These  equations  are  analogous  to  the  generalized  integral  forms  for 
fluid  equations  considering  viscosity  and  heat  conduction.  (3.  5Ala)  says 
there  is  a  mass  source.  In  the  treatment  of  Von  Neumann  and  Hichtmeyer^^ 
this  term  does  not  appear,  since  they  were  particularly  careful  to  maintain  a 
close  analogy  with  the  physical  description.  It  is  one  of  the  aims  of  the 
present  studies  to  investigate  the  present  scheme  with  and  without  this  effect. 

A  first  derivative  approximation  for  at  a  boundary  point  is  required  in 
terms  of  computed  values  'of  U  at  the  interior  points.  Denote  by  P,  A, 
i  s  I,  . . .  N,  the  piston  point,  adjacent  point,  and  subsequent  interior  points. 
Expanding  U  at  P  (J)  with  increments  of  OC  Lr  and  (1  +  )Ar  one  obtains  on 
eliminating  second  order  terms 


1/a  •*'ZOUp 
+'X)  bX" 


(3.5A2) 


Let  this  derivative  expression  be  denoted  by  U.  An  iteration  procedure 
will  be  described  to  evaluate  U  at  the  piston  boundary.  One  obtains  as  a 
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first  approximation  for  p^  in  dimensionless  form 


(3.5A3) 


remembering  that  rarefactions  behind  the  piston  are  to  be  neglected  until  a 
later  date.  This  approximation  uses  the  value  of  r  at  the  previous  time 
(j-1)  and  will  be  improved  once  r^  (j)  is  computed. 

The  remaining  variables  to  be  determined  are  p^r  (t),  tj,  m{r(t),  , 


{r(t).  t  } 


and  r(t).  Two  prdcedures  were  followed  at  this  point. 

As  an  initial  approximation  it  is  assumed  that  the  specific  entropy  is  a 
constant  at  the  piston  (as  it  would  be  in  absence  of  X.  ).  Therefore,  knowing 


Pp  P  p  computed  from 


fp  c:  exp 


(3.  5A4) 


Since  the  density  is  now  determined  over  the  entire  mesh  at  time,  t  , 
one  can  now  compute  the  total  mass.  From  the  equation  of  mass  conservation 
(3.  5Ala)  the  position  of  r  (jAt)  can  be  determined.  However,  there  are  certain 
possible  positions  of  the  piston  for  a  given  At. 

The  position  of  the  piston 


Figure  7 


at  the  time  (jAt)  might  lie  between  the  points  p,  A  of  time  (j-l)At  or  between 
(A,  1),  or  (1,  2),  etc. ,  as  indicated  in  figure  7. 
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Since  a  fixed  meah  is  used  for  all  time,  and  the  quantities  for  all  in¬ 
terior  and  adjacent  points  are  computed  independently  of  the  piston  position, 
it  is  possible  to  integrate  over  a  region  with  too  many  points.  Assume  the 
first  case  described,  use  a  trapezoidal  approximation  for  the  integrals  in 
equation  (3.  SAla)  and  solve  for  cC Lr  to  obtain  the  following  expression  for 
+  1: 


= 


^  j*o  L. 


(3.5A5) 


+.  Wi  - ( I  jj'  -^  ( 1  2oc^")  ?f  j 


IV/-1 


^  2- 


l-ix 


where  the  initial  approximation  oC  ^  at  oCj  ;  is  used. 

Since  it  is  necessary  to  evaluate  V’j  p  at  (j+1),  an  iteration  procedure 
for  was  introduced.  If  0  <  ^  1  then  the  points  over  which  the 

integral  is  evaluated  is  valid  for  a  first  approximation.  Knowing  oC^  +  -^^  and 
the  number  of  points  one  easily  computes  r(j).  It  remains  to  co  pute  m^  (j). 

Using  a  trapesoidal  approximation  once  more  for  equation  (35 Alb),  one 
obtains 


mp '  =  _ 


r 


roil.rrx'r  2oc.v0mr 

j^o  ^  L  /ir 


(3.  5A6) 


^  L,  '=^;s(»  +-'=^:  )Ar 

J  —V. 


1 1,  -+  y  ^  otj.,  hr 
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where  rripj  ^  ^  is  initially  computed  as  r(j+ 1  )-r(j)/ A^t, 

At  this  stage  can  be  computed  from  the  equation  of  state  knowing  p, 
m^,  p  p.  However,  since  the  interior  points  have  computed  with  the 
quasi-viscosity  term  included,  it  would  seem  to  maintain  the  consistency  of 
the  technique  that  E^  be  also  computed  from  an  integral  equation. 

Similarly,  one  obtains 


r 


E 


J--0  L  (y 


(3.5A7) 


?i  oc,(i.<,,Ur-  ]  ' 


N-I 


+ 2E,R.  -  ^  (e-'  7  -  e:;;’ 


^  ot- 


where  Ep'^^  ^  is  initially  computed  from  the  equation  of  state. 

This  determines  all  the  values  for  the  boundary  point.  Equating  the 
equations  (3.  5A3)  and  the  equation  of  state,  one  obtains  the  following  formula 
for  r(t)  at  time  j+ 1 


0.5A8) 
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The  iteration  proceeds  by  repeating  the  sequence  of  calculations  for  p^, 
p  p,  Mp,  Ep,  rp  using  the  latest  computed  values,  averaged  with  the  pre¬ 
ceding  values.  It  is  also  possible  to  avoid  the  use  of  the  assumption  of  con¬ 
stant  specific  entropy.  Again  use  r(j-l)  for  evaluating  Pp.  Now,  only  the 
assumption  of  constant  entropy  at  the  piston,  as  a  first  guess  for  Pp^,  is  used. 


3.5B  Indirect  Methods  for  Boundary  Computations. 

The  alternate  procedure  that  was  tried  was  the  indirect  approach.  This 
method  has  the  advantage  of  simplified  boundary  conditions,  but  requires  a 
compensating  amount  of  computing  time  for  iteration  on  finding  the  true 
boundary, 

R{t)  is  prescribed  for  t  <t<T,  and  assuming  the  piston  path  is  the' par¬ 
ticle  path  of  the  point  initially  at  the  piston,  one  obtains  Up{t)  as  well.  Know¬ 
ing  the  variables  at  the  adjacent  point  A  and  the  piston  point  P,  at  some  time 
t,  the  values  of  the  variables  are  computed  by  linear  interpolation  at  some 
intermediate  point  C  determined  by  the  position  R(t-i-  At) 
figure  8  as 


An.  -  rptt+At)} 


where  U  is  the  vector  p,  m,  E  and  Ar  is  PA.  Using  the  unequally  spaced 
difference  approximations  (3.  4.2)  and  (3.  4.3),  end  applying  them  at  the  point  C 
gives  rise  to  the  following  system 
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ur = -\7.  t  Viivi  x, + Uc" 


(3.  5BZ) 


where 


and 


is  the  vector 


F(U ' ) 


m 


(mjf 


m: 


and  for  the  2nd  component  of 


the  prescribed  value  of  is  used  to  compute  VTl 


j+i 


m 


t 


Having  computed  ^  •)  E 

the  ideal  gas  law. 


X' 
-  fr' 


the  value  is  computed  again  from 


Since  the  equations  for  unequally  spaced  intervals  may  blow  up  as 


o 


j+i 


a  Lagrangean  interpolation  is  performed  to  compute  the 
adjacent  point  when  <  £  Ap  .  This  has  the  tendency  to  greatly 

stabilize  the  unequally  spaced  difference  approximations.  This  feature  was 
not  incorporated  in  the  direct  approach,  though  it  is  now  being  introduced. 

In  the  case,  as  shown  in  figure  9  where  R(t)  may  pass  through  several 
points  for  one  time  step,  two  features  have  been  adopted.  If  only  one  point 
is  passed,  the  newly  computed  adjacent  point  is  omitted.  If  more  than  one 
point  is  passed,  a  refinement  of  the  mesh  spacing  is  performed  by  modifying 


At. 

This  explicitly  defined  numerical  scheme  can  now  continue  to  the  next 
time  step  and  so  on  to  T.  Since  the  magnetic  pressure  has  been  defined  for 
all  t,  comparing  the  computed  pressure  on  R(t)  with  the  prescribed  pressure 
enables  one  to  fit  the  solution  by  modifying  R(t). 

Using  the  results  obtained  from  the  snowplow  model  as  the  first  approxi¬ 
mation  for  R(t),  and  then  interpolating  at  certain  intervals  to  obtain  the 
necessary  agreements  on  the  prescribed  pressure,  leads  to  improved 
approximations  for  the  desired  solutions. 

A  much  more  exact  procedure  for  determining  the  initial  piston  motion 
when  the  piston  is  driven  by  the  magnetic  pressure  is  by  using  the  method  of 
characteristics  prior  to  formation  of  a  shock  wave. 


73 


TN-61-29 


When  this  is  done  for  a  linear  current  relationship  (since  the  shock  form 
is  an  extremely  short  period  of  time)  it  is  found  that  the  piston  velocity  varies 


as  t^;  u  (t)  s  at  where  a  -  (2c  /v)  /  16  P 

p  _  o'  '  o  '  o  o 


The  shock  first 


forms  at  a  time 


ime  t  s  /  2c 

■ 


f  (y)  at  a  distance  from  the  outside  radius 


X.  s  2c 

s  o 


g  (  y  ).  From  characteristic  theory,  therefore,  the 


distribution  of  velocity,  pressure,  and  density  in  the  fluid  between  the  piston 
and  the  point  of  formation  of  the  shock  can  be  readily  derived.  This  enables 
one  to  have  a  more  precise  numerical  procedure  around  small  times. 

Figures  10,  11,  12  show  some  of  the  numerical  results  that  have  been 
obtained  using  these  approximation  schemes. 
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velocity  distribution  fOR.  CONSTANT  RSTON  VELOCITY 

Figure  10 


r,nielssrS 


Ko 


VELOCITY  DiST'n 
FOR  PRE530RE 
.SOURCE  oc  I 
(\-  0005) 


Figure  11 
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Velocity 
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Figure  12 
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EARLY  SOLAR  EVOLUTION 

Or.  Arthur  N.  Cox 
Lo8  Alamos  Scientific  Laboratory 

For  a  number  of  years  several  of  us  at  Los  Alamos  Scientific  Laboratory 
have  been  interested  in  applying  to  astrophysical  problems  some  of  the  tech¬ 
niques  developed  for  bomb  design.  Conditions  in  stellar  interiors  are  some¬ 
what  similar  to  those  in  bombs,  and  in  some  cases  it  has  proved  possible  to 
modify  the  time-dependent  energy  and  mass  flow  equations  so  that  they  can  be 
applied  to  stellar  structure  problems.  With  the  current  ban  on  nuclear  bomb 
testing,  we  now  have  effort  available  to  spend  on  these  astrophysical  applica¬ 
tions. 

There  is  much  work  to  do  in  the  field  of  stellar  structure  and  stellar 
evolution;  and  one  of  the  current  interesting  problems  is  an  understanding  of 
the  complete  history  of  the  sun.  Our  calculations  reported  here  start  at  a 
time  when  probably  the  solar  system  was  condensing  out  of  material  shed 
from  a  collapsing  cloud  of  matter.  We  watch  this  cloud  increase  its  surface 
temperature  from  Just  a  few  hundred  degrees  to  the  present  solar  effective 
temperature  close  to  6,  OOO^K.  Our  future  work  is  to  be  directed  to  the  late 
evolution  of  the  sun  when  the  increased  solar  radius  and  luminosity  will  make 
planetary  conditions  much  different  than  they  are  now.  The  results  I  have  today 
ignore  the  presence  of  the  planets  and  the  rotation  of  the  sun.  Actually  our  data 
apply  only  to  the  rapid  evolution  during  the  first  few  hundred  million  years  of 
the  solar  life. 

To  study  stellar  structure  it  is  necessary  to  consider  a  large  number  of 
physical  processes.  The  equations  that  describe  these  phenomena  are  of  three 
types.  Those  of  the  first  type  in  figure  1  give  the  flux  of  energy  due  to  radia¬ 
tion,  convection,  and  conduction.  The  radiation  flux  is  that  given  when  the 
photons  diffuse  from  one  region  to  another  region  with  only  a  slightly  different 
temperature.  The  Rosseland  mean  absorption  coefficient  or  opacity  appears 
in  this  radiative  flux  expression.  The  other  quantities  are  a,  the  radiation 
constant;  c,  the  velocity  of  light;  p,  the  density;  and  T,  the  temperature  which 
changes  in  our  one-dimensional  calculations  only  in  the  direction  r. 

For  convective  energy  transport  we  use  the  Prandtl  mixing  length 
theory.  The  flux  is  proportional  to  the  1.  5  power  of  the  excess  of  the 
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existing  temperature  gradient  over  the  temperature  gradient  followed  by  an 
adiabatically  expanding  turbulent  element.  The  element  persists  only  for  one 
mixing  length  L,  and  the  mixing  length  is  taken  to  be  1.  5  times  the  local  scale 
height.  In  the  flux  expression  c^  is  the  constant  pressure  specific  heat;  v 
the  turbulent  element  velocity;  P  the  specific  heat  ratio;  constant 

volume  specific  heat;  and  g  the  local  acceleration  due  to  gravity. 

The  conduction  equation  is  the  usual  one  with  v  the  conductivity.  We  note 
that  this  equation  can  be  combined  with  the  radiative  flux  equation  if  a  con¬ 
ductive  opacity  is  properly  defined. 

The  second  set  of  equations  given  in  figure  2  contains  the  hydrodynamic 

equations.  Our  calculations  are  done  in  a  Lagrangian  coordinate  system 

where  we  follow  the  motions  of  mass  shells,  and  it  is  not  necessary  to  write 

an  equation  for  the  conservation  of  mass.  The  conservation  of  momentum 

law  gives  a  force  equation  which  indicates  the  magnitude  and  direction  of  the 

acceleration,  ¥,  of  an  interface  between  two  mass  shells  when  the  pressure 

dP 

gradient,  and  local  gravity  are  not  balanced.  The  Navier-Stokes  hydro¬ 

dynamics  equations  also  include  the  conservation  of  energy  equation  where  Q 
is  the  heat  put  into  a  shell  by  the  flux  equations  or  by  energy  production. 

The  equations  just  described  which  govern  the  flow  of  energy  and  matter 
contain  a  large  number  of  variables,  but  through  various  relations  one  can 
reduce  the  independent  variables  for  a  mass  point  of  fixed  composition  to 
only  the  temperature  and  density.  These  required  relations  are  the  pressure 
and  energy  equations  of  state,  the  opacity  law,  and  energy  production  formulas 
all  of  which  depend  on  the  composition  of  the  matter,  the  temperature,  and 
the  density.  These  relations  are  given  in  figure  3.  Here  b'  is  the  gas  con¬ 
stant  which  includes  electron  degeneracy,  and  and  ^11*6  the  ionization 

and  excitation  energies.  Radiation  pressure  and  energy  have  been  included. 

To  illustrate  some  of  these  relations,  there  are  shown  in  figure  4  pressure 
data  for  a  mixture  of  0.  744  by  mass  of  H,  .  236  by  mass  of  He,  and  .  02  by 
mass  of  heavier  elements  distributed  as  given  on  the  diagrams.  If  the 
electrons  are  nondegenerate  (nuclei  are  always  nondegenerate)  then 

b'  =  b  =  5- 

and  the  pressure  equation  of  state  is  obtained  only  from  the  state  of  ionization 
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of  the  mixture.  At  low  densities  the  b  has  a  plateau  around  a  kT  of  a  few 
volts  because  of  the  complete  ionization  of  hydrogen  and  Hel.  At  higher 
temperatures  the  helium  and  other  rarer  elements  are  ionized  to  produce 
more  particles  for  pressure. 

In  figure  5  the  ionization  energy  of  this  mixture  for  various  temperatures 
and  densities  is  given.  Again,  steps  occur  in  the  variation  with  temperature 
with  constant  density,  because  of  the  successive  ionization  of  the  important 
elements  in  the  mixture. 

Finally  in  figure  6  we  show  the  opacity  of  the  material  at  various  tem¬ 
peratures  and  densities.  A  much  more  detailed  table  of  opacities  than  shown 
has  actually  been  used,  and  the  table  has  been  inserted  in  the  machine  during 
the  computation.  Enough  points  are  included  so  that  linear  interpolation  for 
log  K  versus  log  T  and  log  p  is  accurate  to  10  percent.  Simple  fits  have 
never  proved  very  accurate  over  reasonable  ranges  of  temperature  and  den¬ 
sity. 

The  total  solar  mass  is  divided  into  30  spherical  shells  of  approximately 
equal  thickness  but  not  of  equal  mass.  The  thickness  of  these  shells  now  is 
allowed  to  vary  as  the  usual  time  integration  of  the  energy  and  matter  flow 
requires.  Independent  variables  for  a  shell  of  fixed  mass  are  composition, 
temperature,  and  density.  The  boundaries  between  the  shells  will  take 
positions  such  that  the  net  force  at  an  interface  is  zero. 

From  a  starting  array  of  mass,  composition,  temperature,  density, 
position,  etc.  for  each  mass  shell,  the  flux  equations  and  the  force  equation 
can  be  evaluated  to  give  net  energy  flow  in  the  time  step  /^t  and  the  net 
motion  in  time  At.  There  is  also  to  be  considered  the  thermonuclear  energy 
production  in  the  time  step.  The  energy  equation  of  state  and  the  conservation 
of  energy  equation  can  be  used  to  evaluate  the  temperature  changes  in  the  in¬ 
dividual  mass  shells.  Likewise  the  new  positions  of  the  interfaces  determine 
the  new  densities  in  the  shells.  Any  composition  change  can  also  be  noted, 
and  in  some  cases  this  may  change  the  opacity  law  and  equation  of  state  some¬ 
what.  The  variable  now  can  all  be  advanced  to  be  ready  for  the  next  time 
step. 

After  some  difficulty  with  the  hydrodynamics  as  described  above,  we 
have  converted  to  a  method  which  places  the  interfaces  in  hydrostatic  balance 
at  the  end  of  each  time  step.  Thus  when  going  from  time  n  to  time  n  4-  1  with 
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a  step  of  At, 


i+  1 

+  I 

j=i-l 


i+  1/2 

I 

j=i  -  1 


AT. 

J 


0 


where  for  interface  i  is  obtained  from  the  conservation  of  momentum 
1 

equation,  and  the  partial  derivatives  are  rather  elaborate  analytic  expressions. 
Solution  for  the  Ar  and  AT  values  are  made  by  an  implicit  method  coupled  with 
the  usual  implicit  method  of  computing  radiation  flow.  Never  do  motions  in 
stellar  interiors  get  so  rapid  that  this  assumption  of  constant  perfect  balance 
is  invalid. 

There  is  another  violation  of  the  fundamental  equations  that  we  make,  and 
that  is  in  the  equation  for  the  convective  flux.  The  convection  is  so  efficient 
that  a  very  slight  excess  of  the  existing  temperature  gradient  over  the  adi¬ 
abatic  gradient  can  carry  all  the  energy  required.  Spurious  fluctuations  due 
to  the  integration  procedure  can  then  cause  very  large  excursions  in  the  con¬ 
vective  flux.  Thus,  we  divide  the  convective  flux  by  a  large  factor  so  that  to 
carry  the  flux  required  of  it,  convection  makes  the  existing  temperature 
gradient  somewhat  steeper  than  the  adiabat.  This  steeper  temperature  gra¬ 
dient  only  slightly  affects  the  structure  of  the  star.  There  is  more  room  for 
small  integration  errors  with  this  steeper  gradient,  and  wild  fluctuations  are 
largely  eliminated. 

The  results  I  show  today  are  all  in  the  May  1961  issue  of  Sky  and  Tele¬ 
scope.  For  the  details  one  should  consult  this  article. 

By  manipulation  of  methods  previously  developed  and  by  the  addition  of  a 
few  new  features,  we  have  a  machine  code  capable  of  doing  this  early  evolu¬ 
tion  of  stars  with  both  more  and  less  mass  than  the  sun.  The  later  evolution 
of  stars  is  another  current  problem  which  is  also  being  studied  by  these  time 
dependent  methods. 
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THE  AERONUTRONIC  HOP  PROGRAM  FOR  FLUID  FLOW 

Dr.  R,  A,  Grandy 

Aeronutronic  Division,  Ford  Motor  Company 
The  differential  equations  of  hydrodynamics  are 


dt  p 


dt 


dr  — 

5r  =  *1 


7  P  =  0 

(The  momentum  equation) 

(1) 

^•.”q  =  0 

(The  continuity  equation) 

(2) 

o 

M 

(The  energy  equation) 

(3) 

(The  transformation 
equation) 

(4) 

.  s)  «  g(p,  e) 

(The  equation  of  state) 

(5) 

In  the  above  q  is  the  flow  velocity,  p  the  density,  P  the  hydrodynamic 
pressure,  V  =  1/p  the  specific  volume,  Sthe  entropy,  e  the  energy  density 
per  unit  mass,  and?  =  ^  +  jy  +  is  the  Eulerian  radius  vector  defining  the 
location  of  an  element  of  the  fluid.  We  may  also  write*^  =  iu  +  ^  +  '^w. 

Equations  (1)  to  (5)  are  called  the  Lagrange  equations  of  hydrodynamics. 

In  this  scheme,  y,  and  z  are  considered  functions  of  the  fluid  element  in 
question. 

Finite  difference  methods  for  the  hydro  equations  break  both  the  space  and 
time  continuum  into  finite  intervals.  ^  The  hydrodynamic  variables  are 
sampled  at  various  points  in  space  and  at  various  times.  Let  us  restrict 
attention  to  one-dimensional  flows  for  purposes  of  illustration.  In  this  case 
a  one-  dimentional  spatial  mesh  is  sufficient  to  define  the  problem.  The 
above  equations,  for  example,  are  most  commonly  differenced  in  the  form: 

+  1/2^  (momentum)  (i) 

J  J  (  pA  x) 
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j-1/2  ’■  ^^j-1  ^  J  ^j-l/2[n’‘ 


(continuity)  (2') 


n+  1 

'  J-l/2 


3(!'+*  =  30  <■  U 
j  J  J 


”  _  1  (p*'  + 

j-l/2  r  J-l/2  ^ 

n+  1/2 


pP+  1  )  /v“^  -  V“  ) 

j-l/2'  '  j-l/2  j-l/2' 


(energy)  (3') 
(transformation)  (4') 


_n+  1  / „ n+  1  ^  n+ 1  . 

j-l/2  "  ®^Pj-l/2-  ®j-l/2> 


(equation  of  state) 


(5') 


The  superscripts  n  and  subscripts  j  represent  the  temporal  and  spatial 
mesh  points  under  consideration  respectively.  Time  centering  of  the  above 
equations  has  been  achieved  by  defining  the  velocity  and  coordinate  of  a  mesh 
point  a  half  cycle  apart  in  time.  The  superscript  «•  equals  0,  1,  or  2  for 
slab,  cylindrical,  or  spherical  symmetry.  These  difference  equations  are 
well  known  and,  with  minor  modifications,  are  used  by  the  AFSWC’  in  their 
SHARP  program  and  by  Aeronutronic  in  our  HOP  code.  The  equations  used 
include  an  additional  term  in  the  pressure,  the  von  Neumann-Richtmyer  Q, 
to  allow  the  equations  to  treat  shocks,  but  the  presence  of  this  term  un¬ 
necessarily  complicates  an  illustrative  discussion. 

The  conditional  stability  of  finite  difference  approximations  to  linear 
partial  differential  equations  is  well  known,  having  first  been  discussed  by 
Courant,  Friedrichs,  and  Lewy  in  1928.  ^  A  mathematically  rigorous 
analysis  of  the  stability  of  nonlinear  equations  has  never  been  realized, 
although  nonrigorous  stability  criteria  have  been  established  which  work  in 
practice.  It  is  found  that  as  Ax,  the  mesh  size,  and  At,  the  time  interval, 
separately  approach  zero  in  the  finite  difference  equations,  the  solution  of 
the  difference  equations  does  not  necessarily  uniformly  approach  the  solution 
of  the  differential  equations.  In  most  cases  there  must  be  a  functional  re¬ 
lation  maintained  between  Ax  and  At  as  they  approach  zero.  Normally,  it 
is  found  that  At  must  be  less  than  some  function  of  Ax:  At  <  Ax), 

where  ^  .  represents  parameters  appearing  in  the  differential  equations. 


89 


TN-61-29 


In  our  example,  the  condition  is  that  At  <  Ax  .  where  =  -4^1  is  the 

"  c  dP  /s 

local  sound  speed. 

The  normal  procedure  in  solving  the  finite  difference  equations  is  to  com¬ 
pute  the  quantity  i{»  Ax)  for  every  zone  of  the  mesh,  and,  at  the  end  of  a 
time  cycle,  to  use  the  smallest  of  these  quantities  as  a  time  interval  for  the 
next  time  cycle.  It  is  quite  common  for  f  to  vary  widely  over  the  mesh,  either 
because  cc^  varies  or  because  of  Ax.  In  our  case,  the  sound  speed  might  be 
much  higher  in  some  part  of  the  mesh,  or  the  spatial  interval  might  be 
markedly  decreased  to  give  better  spatial  resolution  in  a  particular  region. 

It  is  normally  necessary,  in  these  cases,  to  calculate  in  the  main  part  of  the 
mesh  with  a  At  much  smaller  than  required  locally  for  either  stability  or 
accuracy,  and  this  can  greatly  increase  the  computing  time  required  for  a 
problem.  Several  distinct  approaches  can  be  made  to  this  problem:  (1)  The 
stability  condition  can  be  relaxed.  A  different  difference  scheme  permitting 
a  larger  At  or  perhaps  even  Ax  can  be  used.  (2)  The  spatial  mesh  can  be 
varied  during  the  course  of  the  calculation  to  minimize  fluctuations  in  the 
time  interval.  (3)  The  local  time  interval  can  be  used  in  advancing  the  spatial 
mesh. 

In  the  case  of  linear  equations,  unconditionally  stable  finite  difference 
schemes  can  be  found  for  most  differential  equations  if  implicit  schemes  are 
considered.  Comparatively  simple  implicit  schemes  can  also  be  obtained  for 
both  Euler  and  Lagrage  hydrodynamic  schemes  which  are  unconditionally 
stable  according  to  the  usual  first-order  stability  analyses.  Unconditionally 
stable  schemes,  however,  do  not  remove  all  restrictions  on  At,  the  time 
interval.  They  allow  it  to  be  chosen  solely  on  the  basis  of  accuracy.  Non¬ 
linear  equations  like  those  of  hydrodynamics  often  permit  the  development  of 
discontinuities:  for  example,  hydrodynamic  shocks.  Contact  discontinuities 
can  also  occur,  even  in  linear  equations.  Normally,  the  difference  equations 
and  spatial  mesh  are  chosen  so  that  these  discontinuities  are  spread  over  two 
or  three  zones,  or  a  real  distance  which  is  as  large  as  can  be  tolerated  for 
the  definition  desired  in  the  problem.  Let  us  follow  Richtmyer's  terminology^ 
and  write  our  differential  equation  as  ^  =  Au.  Most  difference  schemes  are 
based  on  integrating  this  in  the  form 

uH+l  ^  yii  ^  ^(Au)  dt. 

J  J  ^ 
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In  a  simple,  explicit,  forward-time-centered  scheme,  the  integral  is  approx¬ 
imated  by  (Au)^  At,  where  we  have  said  nothing  about  the  space  centering.  An 

J  n't*  1 

implicit  scheme  normally  approximates  the  integral  by  (Au) .  At  (backward 

tn  n*t"  1  ^  ^ 

(Au)j  -I-  J  (centered  time  differ¬ 

encing),  and,  if  the  spatial  differencing  is  done  properly,  unconditionally 
stable  difference  equations  will  usually  result.  In  most  of  our  problems,  Au 
represents  spatial  derivatives  of  u.  As  discontinuities  travel  through  the  mesh, 
these  derivatives,  taken  at  a  point  in  the  mesh,  run  from  infinitesimal  value 
through  finite  value  back  to  infinitesimal  value  as  the  variable  goes  from  its 
initial  to  final  value.  That  is,  (Au)?=  (Au) .  =  0,  We  may  plot  (Au).  as  a 


function  of  time: 


(Au). 


J 


J 


T  is  on  the  order  of  two  or  three  times  Ax/s,  where  s  is  the  propagation 
velocity  of  the  discontinuity. 

If  the  time  interval  is  large  compared  to  t  we  can  easily  make  a  large 
error  in  the  integration  with  any  reasonable  difference  scheme.  If  t”  and 
t”^  ^  are  as  shown,  any  of  the  three  schemes  mentioned  will  give  a  zero  value 
for  the  integral.  For  any  degree  of  accuracy  at  all,  the  time  interval  must  be 
less  than  t  ,  i,  e.  about  1  ( ). 

In  many  cases,  simpler  explicit  conditionally  stable  difference  schemes 
will  permit  a  comparable  At,  so  that  the  more  complex  implicit  schemes  will 
not  necessarily  permit  an  appreciable  reduction  in  computing  time.  Generally, 
the  time  interval  in  an  implicit  scheme  might  be  expected  to  be  determined  by 
that  region  of  the  problem  where  quantities  are  rapidly  varying,  for  example, 
in  the  vicinity  of  a  shock  front.  The  requirement  may  be  that  At  .^h(  j, 

Ax)j  where  j  is  a  point  in  the  shock  region.  An  implicit  scheme  can  be  very 
valuable  in  cases  where  h  has  a  reasonable  value  for  most  of  the  problem. 

That  is,  even  if  a  small  portion  of  the  mesh  is  finely  zoned  for  better  detail, 
this  portion  of  the  mesh  will  only  affect  the  time  interval  when  a  shock  passes 
through  it.  Otherwise,  it  is  not  artificially  slowing  the  computing  merely  to 
keep  stability  in  its  vicinity.  Another  approach  which  has  been  exploited  in 
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solving  other  types  of  initial  value  problems,  like  diffusion  problems,  is  the 
use  of  higher  order  schemes  in  the  spatial  differencing  for  the  purpose  of  in¬ 
creasing  the  spatial  interval  necessary  for  a  given  desired  definition.  Such 
schemes  have  been  used  occasionally  in  hydrodynamics  to  reduce  truncation 
errors  in  cases  where  the  mesh  size  changes  rapidly  over  the  mesh,  but  have 
not  been  seriously  exploited  in  the  above  sense. 

There  are  many  problems  of  interest  where  an  explicit,  conditionally 
stable  scheme  will  have  a  local  time  interval  which  varies  widely  over  the 
mesh  and  where  the  time  interval  required  for  accuracy  in  an  unconditionally 
stable  scheme  will  always  be  smaller  than  desired.  In  purely  hydrodynamic 
problems,  it  has  been  found  to  be  very  profitable  to  use  the  local  time  interval 
in  advancing  the  spatial  mesh.  The  following  logical  scheme  permits  each 
zone  of  the  mesh  to  have  its  own  time  interval,  and  advances  the  time  by  an 
arbitrary  amount  per  machine  cycle. 

At  the  beginning  of  a  cycle,  all  points  of  the  mesh  exist  at  the  same  time. 
Generally,  after  a  further  sweep  through  the  mesh,  the  points  will  exist  at 
different  times,  A  point  is  not  accelerated  if  any  of  its  neighbors  have  not 
yet  been  advanced  as  far  in  time  as  the  point  in  question.  If,  however,  the 
point  is  at  a  time  less  than  or  equal  to  that  of  all  neighbors,  it  is  accelerated. 
The  neighbors  are  relocated  temporally  by  interpolation  to  the  time  of  the 
point  in  question.  Hydrodynamic  quantities  are  calculated  at  the  time  of  the 
point  and  the  point  is  moved  using  as  a  time  interval  the  smallest  of  the  neigh¬ 
boring  intervals. 

Various  special  cases  may  arise.  If  any  of  the  neighbors  is  to  be  ac¬ 
celerated  at  the  same  time  as  the  point,  the  interpolation  routine  is  redundant 
and  is  bypassed  for  that  point.  If  the  time  interval  is  such  as  to  advance  the 
point  beyond  the  end  of  cycle  time,  it  is  adjusted  to  just  reach  this  time.  No 
zone  is  allowed  a  time  interval  greater  than  the  cycle  time  interval.  This 
allows  inactive  points  or  points  With  a  large  natural  time  interval  to  be  ad¬ 
vanced  to  the  end  of  cycle  time  in  one  sweep.  If  any  of  the  neighbors(of  a 
point  which  has  been  advanced  to  the  end  of  cycle  time  with  the  cycle  time 
interval)  have  not  yet  been  advanced  to  this  time,  the  inactive  point  is  backed 
up  to  the  earliest  such  time  and  reaccelerated.  This  test  automatically 
activates  points  when  a  disturbance  reaches  them. 
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The  above  scheme  has  been  applied  to  the  Lagrange  difference  equations 
for  one-  and  two-dimensional  flow.  We  shall  briefly  describe  the  Aeronutronic 
HOP  program  for  one-dimensional  flow.  This  program  was  originally  de¬ 
veloped  for  use  on  Bureau  of  Ordnance  contract  NOrd  17945,  concerning  the 
detonability  of  propellants. 

The  HOP  program  is  best  explained  by  reference  to  the  flow  chart. 

(Figure  1.  )  Location  I  is  the  entry  point  to  begin  a  macro  cycle.  Part  of  the 
left  hand  boundary  conditions  are  set,  registers  are  initialized,  and  the  pro¬ 
gram  proceeds  to  A,  At  A,  the  time  at  which  the  coordinates  of  point  j  exist 
is  compared  with  the  time  for  point  j  +  1.  If  t^  >  t^  ^  the  program  exits 
to  make  an  activity  test  at  location  M.  If  t^  <  ,  it  proceeds  to  acceler¬ 

ate  the  point  unless  it  is  already  at  the  end  of  cycle  time,  in  which  case 
another  activity  test  is  made.  If  L  <  ^  x^^  is  recomputed  at  the  time 

of  point  j  under  the  assumption  that  j+  1  was  accelerated  with  a  constant 
acceleration.  If  t^  s  fbis  interpolation  is  bypassed  and  we  proceed 

directly  to  calculate  hydrodynamic  quantities  for  zone  j+Vz  at  the  time  of  the 
point  j.  If  tj  s  tj^  and  if  it  turns  out  that  we  may  also  accelerate  point 
j+  1  during  the  same  sweep,  we  will  not  need  to  calculate  hydrodynamic 
quantities  for  the  left  hand  zone  when  doing  j+1.  Sense  light  1  is  turned  on  to 
inform  the  program  that  this  is  so. 

Location  C  is  the  exit  from  calculating  hydrodynamic  data  for  zone  j+i/fe. 

At  this  point,  sense  light  2  is  checked.  If  it  is  on,  it  means  that  hydrodynamic 

data  for  zone  j-i/z  at  time  j  already  exists.  If  it  is  off,  t^  is  compared  with 

tj  J,  and  in  either  case  hydrodynamic  data  are  calculated  for  zone  j-Vz. 

Location  D  is  the  exit  from  the  hydrodynamic  calculations  for  zone  j-V^.  The 

program  transfers  directly  to  D  if  sense  light  2  is  on.  At  D,  the  smallest 

of  the  two  time  intervals  At  . . and  the  cycle  interval  At  is  chosen  as  the 

J±V2  C 

time  interval  for  the  point  j.  The  time  is  advanced  and  checked  against 
N+  1 

t  .  If  the  new  time  is  less  than  the  end  of  cycle  time,  sense  light  3  is 
turned  on  as  an  indication  that  at  least  one  point  has  not  yet  been  advanced  to 
the  end  of  cycle  time.  The  point  is  then  accelerated  and  a  new  velocity  and 
coordinate  are  calculated. 

Sense  light  1  is  on  if  the  old  L  =  lu  this  case,  sense  light  2  is 

turned  on  and  the  program  transfers  to  B  where  it  proceeds  to  process  the 
next  point.  If  sense  light  1  is  off,  t^  was  less  than  and  normally  the  pro- 
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N+  1 

gram  skips  and  proceeds  to  test  point  j  +  2  unless  ^  !  If  this  is  true, 

j  is  stepped  by  one  and  control  passes  to  the  activity  test  M.  Normally  j  is 
stepped  by  two  and  the  transfer  is  to  location  G  where  j  is  checked  against 
JMAX.  If  j  >  JMAX,  the  program  transfers  to  F  to  begin  a  new  sweep  through 

the  mesh.  If  j  £  JMAX,  t.  is  compared  with  t.  ,  .  If  t .  >  t .  . ,  it  is  compared 

N+1  ^  N+ 1 

with  t  .  If  these  are  equal,  the  activity  test  is  made.  If  L  <  t  .  j  is 

stepped  by  one  and  the  program  returns  to  G.  If  t^  <  t^  the  point  may  be 
accelerated  depending  on  If  j  •  JMAX,  the  program  transfers  to  H 

where  the  right-hand  boundary  conditions  are  set  and  an  activity  test  made. 

The  activity  test  at  location  M  first  compares  At^  against  At^,  the  cycle 
interval.  If  Atj  <  At^,  the  point  has  actively  reached  the  end  of  cycle  time 
and  the  program  proceeds  to  step  j+ 1  and  does'the  next  point.  If  the  two  in¬ 
tervals  are  equal,  the  program  selects  the  smaller  of  times  t.  ,  and  t-. 

N+ 1  J“  ^  ^ 

If  the  smaller  of  these  is  equal  to  t  ,  the  program  proceeds  to  do  the  next 
point.  Otherwise,  the  point  j  is  relocated  to  a  time  t,  where  t  refers  to  the 
smaller  of  the  adjacent  times  and  At  is  the  corresponding  time  interval.  Point 
j  is  then  accelerated. 

Location  F  begins  a  sweep  through  the  mesh.  If  sense  light  3  is  on,  at 

least  one  point  is  not  yet  at  the  end  of  cycle  time  and  control  passes  through  I. 

N+ 1 

If  sense  light  3  is  off,  the  entire  mesh  has  been  advanced  to  time  t 
Various  end-of-cycle  tests  are  made  and  the  cycle  time  is  stepped.  The  pro¬ 
gram  then  transfers  to  I  to  begin  a  new  macro  cycle. 

One  further  departure  is  made  from  the  standard  difference  equations. 

Since  ^  s  g  the  time  centered  energy  equation  (3')  is  an 

implicit  equation  as  e*'^  ^  occurs  on  the  right  hand  side.  If  P  is  a  linear 
function  of  e,  which  would  be  the  case  for  a  polytropic  gas,  equation  (3')  can 
easily  be  solved  for  e*^^^  and  an  explicit  time  centered  equation  obtained. 

Even  a  second  order  equation  has  some  truncation  error.  In  many  problems, 
the  y^PdV  must  be  performed  over  a  change  in  density  of  several  orders  of 
magnitude,  arriving  at  a  final  pressure  many  orders  of  magnitude  less  than  the 
original  pressure  which  must  be  given  with  an  accuracy  of  at. least  several 
percent  relative  to  its  final  value.  It  is  impossible  to  get  the  desired  accuracy 
using  equation  (3')  ,  For  this  reason,  we  use  the  equation  of  state  in  the  form 
P  .  f  (p,  S).  During  the  adiabatic  expansion,  S  is  constant,  or  at  least 
slowly  varying.  For  a  polytropic  gas,  P  ■  A(S)p^.  This  is  our  starting  point 
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in  integrating  the  energy  equation  and  equation  of  state.  Basically,  with 
corrections  for  the  von  Neumann- Richtmyer  Q  and  for  other  source  or  sink 


terms,  the  pressure  equation  is  written  as 
,  n+ 1  Y 

pn«  .  pn 


where  y  a  slowly  varying  function  of  P  and  p  which  incorporates  the  degree 
of  ionization  of  the  material. 
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SOME  COMMENTS  ON  CONVAIR'S  WORK 
ON  NUMERICAL  METHODS  OF  FLUID  FLOW  PROBLEMS 


W.  F.  Brown 
C.  G.  Davis 
D.  A,  Hamlin 
A.  H.  Schainblatt 

Convair,  San  Diego 


1.  INTRODUCTION. 

A  discussion  of  some  of  the  work  on  Numerical  Methods  of  Fluid  Flow 
Problems  that  has  been  carried  out  by  the  Physics  Section  of  Convair-San 
Diego  will  be  given.  This  discussion  will  include  general  philosophy,  des¬ 
cription  of  programs,  techniques  used,  and  results. 

In  order  to  avoid  any  semantic  difficulties,  let  us  state  that  when  we  talk 
about  one-  or  two-dimensional  hydrodynamics  we  mean  one  space  variable 
and  time  or  two  space  variables  and  time,  respectively.  For  all  other 
problems,  it  will  be  assumed  that  they  are  time  independent  unless  otherwise 
noted;  thus,  when  we  talk  about  a  one -dimensional  non-equilibrium  flow,  we 
mean  a  time  independent  flow  in  one  space  variable.  The  aireas  which  we  will 
cover  are  one -dimensional  and  two-dimensional  hydrodynamics  and  plane  and 
axially  symmetric  non-equilibrium  flows. 

2.  GENERAL. 

It  seems  appropriate  to  give  a  few  general  statements  of  our  philosophy 
regarding  machine  codes.  We  attempt  a)  to  prepare  our  codes  so  that  they 
are  not  the  personal  property  of  the  individual  scientists  who  developed  them; 
b)  to  use,  whenever  possible,  tested  techniques;  and  c)  to  write  our  codes  in 
such  a  way  that  they  may  be  generalized  but  still  specific  enough  to  do  the 
job  at  hand  in  a  reasonable  time.  Furthermore,  whenever  possible,  analyti¬ 
cal  and  numerical  techniques  are  used  in  a  complimentary  way. 

Since  we  will  have  need  of  them  later,  let  us  now  write  the  hydrodynamic 
equations  in  Lagrangian  coordinates.  In  order  to  do  this,  we  now  represent 
the  original  Lagrangian  coordinates  of  a  fluid  element  as  small|xj|.,i  ■1,2, 

3.  and  the  Eulerian  coordinates  as  large|x^|  ,  i  ■  1,  2,  3. 

The  Eulerian  coordinates  which  may  be  considered  as  functions  - 
X^(x^,  t),  i  ■  1«  2,  3,  give  the  positions  at  time  t  of  the  fluid  elements  which 


97 


TXM-61-Z9 


were  originally  at  x^. 

Under  the  above  definitions  the  law  of  conservation  of  mass  becomes 


v.  3  tx,3 


(2.  1) 


where  V  s  1/p  is  the  specific  volume  of  the  fluid  element.  is  the  original 


specific  volume  in  the  Lagrangian  system  and  ^ 
Jacobian  of  the  transformation 

We  may  write  the  equation  of  motion  as 


is  the 


=  -V 


1,2,3 


(2.2) 


Our  basic  set  of  differential  equations  are  completed  by  the  conservation 
of  energy,  which  may  be  written  as 


JE  =-pdv  +  5 


(2.3) 


where  E  is  the  specific  internal  energy  of  the  fluid  and  S  is  a  quantity  which 
contains  sink  as  well  as  source  terms. 

In  order  to  complete  the  above  set  of  equations,  an  equation  of  state  is 
needed.  One  may  use  for  this 


p  =?(fEy) 


(2.  4) 


3.  ONE-DIMENSIONAL  HYDRODYNAMICS. 

By  specializing  the  above  equations  to  one  space  dimension  and  taking 
central  differences,  we  obtain  the  difference  form  of  the  one -dimensional 
hydrodynamic  equations  when  S  >  0. 
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Here  P  is  the  fluid  pressure,  a  is  the  "shock  width  constant,  "  R  ■  X^,  q 
is  the  pseudo-viscosity  pressure  term  and  0Cs  1,  1,2,  for  slab,  cylinder  and 
sphere  respectively.  If  we  now  consider  the  cas«  of  an  ideal  gas,  we  can 
write  (3.  4)  as 


(3.6) 
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Also,  for  our  itability  criteria,  we  have  uaed  the  White  number 


(3.7) 


where  C  is  the  velocity  of  sound. 

Let  us  now  consider  a  spherical  piston  expanding  into  an  ambient  atmos¬ 
phere  and  assume  that  the  flow  field  external  to  the  piston  is  represented  by 
the  one -dimensional  hydrodynamic  equations.  We  will  now  consider  three 
models  for  the  internal  flow  field. 

The  first  will  take  the  piston  as  a  thin  shell  with  constant  mass  and  given 
initial  energy.  Interior  to  the  mass  shell,  it  will  be  assumed  that  Ps  0.  For 
this  case,  the  equations  governing  the  motion  of  the  shell,  which  we  represent 
as  a  mass  point,  are  given  in  difference  form  as 


(3.8) 


R~‘-.  MU.r‘  *  Kl 

whereP^  (the  external  pressure)  is  found  from  the  hydrodynamic  equations 
and  M  is  the  total  mass  of  the  piston. 

In  the  other  two  models  which  we  have  considered  to  date,  an  internal 
structure  is  assumed  for  the  piston.  The  flow  internal  to  the  piston  is  con¬ 
sidered  to  have  a  uniform  initial  density  which  is  sufficiently  low  that  particle- 
particle  collision  can  be  neglected  and  an  expansion  velocity  initially  given  by 

U  = 

In  the  first  case,  it  is  assumed  that  the  interface  is  a  membrane  and  the 
internal  particles  reflect  elastically  from  the  membrane.  Under  these 
assumptions,  the  equation  governing  the  motion  of  the  piston  is  given  by 

(i  -  X-  i  X*)  =  <  1-*  (i  -  K")  ) 

(3.  9) 
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ts  = 


S  ■  V4TrPo 


f  = 


1  - 
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Here  f  and  X  denote  the  fractional  lag  in  velocity  and  position,  respectively, 

of  the  actual  interface  behind  the  unretarded  interface,  and  t  is  a  scale  time. 

s 

The  other  case  considered  assumes  that  when  a  particle  hits  the  interface 
it  sticks  to  it.  In  this  case,  the  governing  equation  becomes 


where  the  variables  are  defined  as  before,  except  that  OC^  is  replaced  by 

Based  on  the  above  equations,  a  series  of  codes  have  been  written  and  a 
number  of  cases  computed. 

The  first  program  solved  the  case  of  a  spherical  piston  expanding  into  a 
constant  ambient  atmosphere.  In  this  code,  a  fixed  number  of  points  were 
maintained  in  the  flow  field  and  only  the  points  which  were  influenced  by  the 
piston  motion  were  calculated.  In  figure  1,  we  see  the  results  of  some  typical 
calculations.  In  this  figure,  the  circled  points  denote  calculations  made  for 
yz  5/3  and  triangles  for  2.  As  is  clearly  seen,  the  variation  of  9* 
has  a  very  small  effect  on  the  solution.  Other  calculations  have  been  made 
for  different  initial  conditions  and  the  results  are  similar  to  those  shown  in 
figure  1.  Typical  pressure  and  velocity  profiles  obtained  from  our  machine 
calculations  are  shown  in  figures  2  and  3  respectively. 
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Figure  1.  Velcx:lty  vs.  Distance  for  Expansion 
of  a  Mach  6  Spherical  Shell  Into  a  Uniform 
Ambient  Atmosphere:  Conq^iarlson  of  y  «  5/3 
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G.  I.  Taylor^  has  obtained  a  solution  for  a  piston  expanding  with  a  constant 
velocity  in  an  ambient  atmosphere.  From  this  solution,  one  can  see  that  over 
a  wide  range  of  velocity  we  can  represent  the  pressure  on  the  piston  in  the  form 

where  c  is  the  velocity  of  sound  in  the  ambient  medium. 

This  can  also  be  written  as 


(3.  12) 

where 

-  r  i^T 

To  determine  F  •  we  make  use  of  numerical  results  obtained  from  the 

expansion  of  a  spherical  piston  with  constant  velocity.  These  results  are 
illustrated  in  figures  4,  5,  6. 

Using  this  form  of  the  pressure,  one  can  obtain  an  analytic  solution  for 
the  motion  of  the  mass  shell  piston.  The  results  are  given  in  figure  7  and  a 
comparison  is  made  with  a  hydrodynamic  calculation  in  figure  8. 

Because  of  the  applicability  of  the  fitted  pressure,  it  was  decided  to  solve 
the  other  two  piston  models  using  this  form  of  the  pressure,  but  analytical 
solutions  cannot  be  obtained  as  in  the  case  of  the  shell  model,  except  for  the 
special  case  of  a  constant  pressure  for  which  a  power  series  solution  can  be 
obtained.  When  the  pressure  equation  (3.  12)  is  substituted  into  equations 
(3.  9)  and  (3.  10),  we  have  ordinary  differential  equations  for  the  interface 
velocity.  A  program  was  prepared  which  integrates  these  equations  and 
results  of  typical  calculations  are  given  in  figure  9.  Also,  as  part  of  these 
codes,  the  partition  of  energy  was  calculated  and  results  of  a  typical  case  are 
shown  in  figure  10,  where  denotes  the  streaming  energy,  is  the  kinetic 
energy  of  the  shell  in  the  sticky  model,  E^  is  the  heat  produced  in  the  shell 
because  of  the  sticking,  and  E^  is  the  work  done  against  the  total  back  pressure; 
E^  and  E^  are  not  applicable  to  the  reflecting  model.  For  the  sticking  model, 
the  total  energy  is  clearly  conserved  to  the  accuracy  of  the  calculation.  The 
total  energy  for  the  reflecting  model  is  very  difficult  to  calculate  and  has  not 
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Figure  4.  Pressure  at  Surface  of  Spherical 
Piston  Expanding  at  Constant  Velocity  for 
Several  Values  of  v  • 
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Figure  4 


Figure  Pressure  at  Surface  of  Si^erlcal 
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Figure  5 


Figure  6.  Ratio  of  Sboek  Radius  to  Plstou 
Radius  for  Qptasrlesl  PlstOD  liqpandliig  at 
Constant  Velocity  Uj^/C  ■  3  for  Several 
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Figure 
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Figure  10.  Distribution  of  Energy  vs.  Badlus  for  Reflecting 
and  Sticking  Models  of  Bicpanslon  of  a  Mscb  7*77  Uniform 
Stphezlcal  doud  into  a  uniform  Ambient  Atmosphere. 
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yet  been  done. 

In  order  to  check  these  results,  a  hydrodynamic  program  has  been  pre¬ 
pared  which  uses  equation  (3.  9)  and(3.  10).  Preliminary  results  from  this 
code  indicate  that  the  pressure  given  by  equation  (3.  11)  is  still  a  good  approx¬ 
imation.  In  figure  11,  a  comparison  is  given  between  all  three  models  for  the 
pressure  given  by  equation  (3.  11). 

Programs  have  also  been  prepared  which  calculate  the  expansion  into  a 
non-uniform  atmosphere  for  the  three  models  under  consideration.  Presently, 
calculations  are  being  made.  It  should  be  noted  that  with  our  present  codes  we 
could  calculate  piston  motion  when  the  internal  flow  is  governed  by  the  hydro- 
dynamic  equations. 

In  this  work,  we  have  an  example  of  analytical  and  numerical  methods 
complementing  each  other;  also,  of  a  machine  program  being  flexible  enough 
to  permit  extensions  and  yet  give  useful  results  at  each  stage  of  development. 

Other  problems  in  this  area  in  which  we  are  actively  engaged  are  hydro- 
dynamic  flow  with  radiation.  Our  present  work  is  limited  to  radiation  diffusion. 
Thus,  the  only  difference  that  occurs  in  our  basic  equations  is  in  the  energy 
equation  which  can  be  written  in  differential  form  as 

= -pdv  +  v-(Kve‘‘) +S 

where  K  is  a  function  of  the  opacity,  6  the  temperature,  and  S  a  source  term. 
Our  present  code  is  so  constructed  that  the  energy  equation  is  solved  in  a 
subroutine,  thus  permitting  different  differencing  techniques  to  be  used  in  the 
diffusion  term  and  use  of  different  methods  of  solving  the  resulting  nonlinear 
equations.  In  the  present  code,  the  Crank-Nickolson  method  is  used.  This  is 
an  implicit  method  which  replaces  the  diffusion  term  at  n  +  1/2  by  the  average 
of  term  at  n  and  n  +  1.  The  radiation  coupling  is  done  by  first  letting  the 
radiation  flow  from  n  to  n  +  1  with  the  hydrodynamics  held  fixed,  and  then  using 
the  temperature  at  n  +  1,  compute  E  and  P.  Now  a  number  of  other  methods  of 
coupling  can  be  used  but  the  method  presently  being  used  has  been  successfully 
tested  and  used  on  similar  problems. 

In  our  radiation-hydrodynamics  effort,  we  have  an  example  of  a  code  which 
has  developed  from  a  fairly  simple  code  to  one  of  some  complexity  in  which  the 
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first  methods  being  used  are  tested  but  the  code  is  so  constructed  that  basic 
changes  can  be  made  in  the  code  without  major  reprograming. 

4.  TWO-DIMENSIONAL  HYDRODYNAMICS. 

In  this  area,  we  have  been  interested  in  both  pure  hydrodynamic  as  well 
as  radiation-hydrodynamic  calculations. 

In  order  to  achieve  our  goals,  we  first  prepared  a  hydrodynamic  code 
modeled  after  the  Los  Alamos  Magee  code.  Thus  the  difference  equations 
used  are 
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AV”*'  =  V"*'  -V" 


(4.8) 

(4.  9) 

(4.  10) 

(4.11) 
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The  acceleration  terms  are  those  obtained  by  the  so-called  Force  Gradient^ 
Method.  This  code  has  the  feature  that  the  form  of  the  acceleration  terms 
can  easily  be  changed  as  well  as  the  equations  of  state;  also,  the  overall  logic 
is  such  as  to  minimize  the  time  of  calculations. 

In  order  to  develop  our  radiation  hydrodynamics  code,  we  first  prepared 
a  radiation  diffusion  code  which  can  be  coupled  to  the  hydrodynamics  in  the 
same  way  as  for  the  one-dimensional  case.  The  basic  equation  solved  in  the 
radiation  code  is  the  energy  which  when  written  in  integral  form  becomes 


(4. 16) 


where  £  *  £  (  d,V  )  is  the  specific  internal  energy,  i.  e.  ,  per  unit  original 
volume,  ([)  is  the  flux  of  radiation  across  the  surface  S  bounding  the  element 

^  o* 

By  differencing  these  equations  by  the  implicit  alternating  direction 
method,  i.  e. ,  for  one  time  step,  the  equations  are  differenced  implicity  in 
the  X  direction  and,  on  the  next  time  step,  implicitly  in  the  Y  direction.  We 


obtain 


A  AO’  "i" 


.n+'/t 


n+Zz 


17) 


=  'DJ' 


X+«jX-*l  T+l,T+l 


(4.18) 


A  r\  T  q 

coefficients  ^  ^rv>jT+\  nre  known  functions 

of  known  quantities  at  time  n. 

These  equations  are  solved  by  introducing  the  well  known  recursion 
relationship 
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-a"”  +  b"" 

^x-i,i+\  r-i,T-v\  x-n,r+i 

By  substituting  these  into  the  above  equations,  one  obtains 

ant^i 
X+ljT+j 
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(4.  19) 
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(4.21) 


and  a  similar  set  for  the  Y  implicit.  Using  these  and  the  boundary  conditions, 
one  is  able  to  solve  the  necessary  equations.  The  code  used  to  solve  these 
contains  the  same  feature  as  our  one -dimensional  code  and  the  coupling  is 
done  in  the  same  way  as  before. 

7 

Here  again,  we  have  used  tested  methods  but  still  allow  for  other  methods 
to  be  easily  used,  and  have  developed  our  code  in  a  step-wise  manner. 

5.  NON-EQUIUBRIUM  FLOWS. 

Let  us  now  briefly  turn  to  the  area  of  non -equilibrium  gas  flows.  We 

4 

have  chosen  as  our  model  that  given  by  Wood  and  Kirkwood  with  the  further 
assumption  that  we  have  local  thermodynamic  equilibrium  between  classes 
of  degrees  of  freedom  as  well  as  within  these  classes.  These  equations 
could  also  have  been  obtained  by  specializing  the  Kirkwood-Crawford^ 
equations  to  an  inviscid  nondiffusing  gas.  The  equations  are 

vH  -  v/^ 


4  VP  =0 
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d*t 


(5.3) 


1  =  -  £  A.<.Fc^  (5.4) 

c5t  “ 


The  density  derivative  has  been  eliminated  from  the  mass  equation  by  use 
of  the  differential  form  of  the  caloric  equation  of  state 


(5.5) 

where  we  have  defined 

a;  =  [S%/c^p2  (5. 6) 


(frozen  velocity  of  sound) 
(frozen  expansion  coefficient) 


(frozen  specific  heat  at 
constant  pressure) 


(5.7) 

(5.8) 

(5.9) 

(5.  10) 


Also  Pg  is  tbe  specific  chemical  potential  of  specie  s  and  the  specific 
stoichiometric  coefficient  for  the  s  specie  in  the  acth  reaction. 

In  order  to  complete  our  set  of  equations,  we  introduce  the  equation  of 
state  for  an  ideal  mixture 


(5. 11) 


Here  Z  is  the  compressibility. 

The  above  equations  govern  the  flow  in  non-equilibrium.  The  equations 
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for  frozen  flow  are  obtained  simply  by  setting  V(x.  •  0.  It  can  be  shown^ 
that  the  equilibrium  equations  can  be  obtained  from  the  frozen  equations  by 
replacing  the  frozen  velocity  of  sound  by  C^,  the  equilibrium  velocity  of 
sound,  and  evaluating  all  of  the  thermodynamic  functions  from  equilibrium 
thermodynamics.  It  is  easily  seen  that  these  basic  equations  are  hyperbolic 
in  character  and  that  the  difference  between  equilibrium  and  non-equilibrium 
equations  is  only  in  the  coefficients  and  added  flux  type  terms.  Thus  one  can 
use  the  usual  tools  for  calculation  of  hyperbolic  equations.  limiting  ourselves 
to  plane  and  axially  symmetrical  flows,  we  find  the  characteristic  relations 
to  be 


"5r  i  ~  L.  c;«  -  U.» 


(S.  12) 


where  the  plus  sign  denotes  the  right-running  characteristic  and  the  minus 
sign  the  left-running  characteristic.  Also,  the  streamlines  are  character¬ 
istic,  i.  e. , 

6^  (5.  IJ) 

The  compatibility  equations  become 

^Co(^aci Y  -  v6ub)  +  -  ^o)  F-  -  \^\/]  (Jp 

(5. 14) 

where  b  s  1  for  axially  symmetrical  flows  and  0  for  plane  flows.  Another 
method  is  to  finite  difference  the  basic  equations  in  such  a  way  that  the  cal¬ 
culation  is  carried  out  using  the  proper  domain  of  dependence  determined 
by  the  characteristics.  This  has  been  done  but  will  not  be  listed  here.  In 
order  to  carry  out  calculations  using  the  above  equations,  one  needs  the 
thermodynamics  which  for  equilibrium  can  be  found  in  essentially  three 
different  ways,  (a)  table  look  up,  (b)  analytical  fits  to  the  tables,  and  (c) 
calculation  of  the  basic  equilibrium  equations.  It  has  been  shown*^  that  for 
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a  wide  range  of  interest  air  can  be  represented  as  a  simple  gas  governed  by 
the  reactions 

Oa  ^  20 
^  2  hJ 

No  ^  0  +  KJ 


Also,  all  the  needed  thermodynamic  quantities  can  be  expressed  in  terms  of 
G,  the  dimensionless  internal  energy  Z  and  all  their  derivatives. 

Using  these  concepts,  subroutines  have  been  prepared  for  calculation  of 
the  thermodynamics  for  equilibrium  and  non- equilibrium.  The  basic  difference 
in  these  subroutines  is  that  for  equilibrium  the  species  concentrations  are 
calculated  from  the  equilibrium  equations  which  are  nonlinear  algebraic  equa¬ 
tions,  while  for  the  other  cases  they  are  obtained  from  the  kinetics  which  are 
solved  as  part  of  the  basic  flow  equations.  Codes  have  been  prepared  which 
use  these  subroutines  and  the  basic  equations  for  calculating  flow  fields.  It 
has  been  determined  that  it  is  wiser  to  develop  separate  programs  for  equili¬ 
brium  and  non -equilibrium  rather  than  to  try  to  calculate  all  cases  with  one 
code. 

The  one  basic  difficulty  introduced  by  non-equilibrium  flows  over  and  above 
determination  of  reaction  rates  is  the  limitation  placed  on  the  mesh  size  by 
the  kinetics.  It  has  been  found  in  some  cases  that  the  mesh  size  must  be 
taken  such  that  it  corresponds  to  a  flow  time  of  a  nanosecond.  Consequently, 
large  amounts  of  machine  time  must  be  expended  for  a  calculation.  At  the 
present,  we  are  investigating  various  methods  of  reducing  the  calculating  time. 
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PROBLEMS  OF  RADIATION  TRANSPORT 
IN  HEATED  AIR 

Dr.  R.  K.  M,  Landshoff 
Dr.  R.  E.  Meyerott 
Lockheed  Missiles  and  Space  Division 

Nuclear  explosions  in  the  air  can  be  strongly  influenced  by  the  emission 
and  absorption  of  radiation  which  transfers  energy  from  one  point  of  the  fire¬ 
ball  to  another  or  to  points  in  the  exterior.  The  flow  of  radiation  can  be 
expressed  in  terms  of  the  intensity  I^,  where  I^  cos  PdOdvdtdA  is  the  energy 
passing  through  an  area  dA,  during  a  time  dt,  in  a  frequency  interval  dv  and 
coming  from  an  element  of  solid  angle  d  Q  around  a  ray  which  makes  an 
angle  P  with  the  normal  to  dA.  In  general,  the  intensity  I^  is  a  function  of 
position  and  time,  direction  of  ray,  and  frequency  or  wavelength. 

Along  a  given  ray,  the  intensity  changes  in  a  manner  which  is  determined 
by  emission  and  absorption  of  the  radiation  by  the  material  through  which  it 
passes.  Such  changes  are  calculated  from  the  equation  of  radiative  transfer. 
If  the  material  is  in  local  thermodynamic  equilibrium,  the  equation  of  trans¬ 
fer  is  ^ 


where 


dl  V 

ds 


K 


(B, 


(1) 


(1  -  exp  (  -  hv/kt)  ) 


is  defined  in  terms  of  the  absorption  coefficient  K  (p,  T);  s  measures  the 


distance  traversed  by  the  radiation  along  the  ray;  ^ 


2h  V 


exp  ( _ )-l 

kt 


The  mathematical  character  of  equation  (1)  is  well  known  and  leads  to 
solutions  where  I^  ,  tends  towards  .  If  the  temperature  remains  constant 
in  a  large  enough  region  of  space  rather  than  changing  from  point  to^point, 
i.  e.  ,  if  one  has  complete  rather  than  local  thermodynamic  equilibrium,  B^ 
will  also  be  constant  and  I  will  become  equal  to  B  .  This  is  the  situation 
within  a  "Hohlraum"  or  blackbody,  and  B^  is  thus  the  blackbody  intensity. 

Equation  (1)  treats  the  radiative  transfer  problem  for  each  instant  as  if 
the  density  and  temperature  distributions  were  held  fixed.  Such  an  adiabatic 
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treatment  is  justified  because  of  the  high  velocity  of  light. 

If  the  absorption  coefficient  ^  1*  known  as  a  function  of  p  and  T,  one  can 
formally  solve  equation  (I)  in  terms  of  the  instantaneous  spatial  distributions 
of  p  and  T.  To  write  this  solution  down,  it  is  convenient  to  introduce  the 
optical  depth  of  a  general  point  P  on  the  ray. 


Tv 


i  “v'*' 


The  upper  limit  is  a  point  beyond  the  edge 
where  the  ray  emerges  from  the  absorbing 
region  into  the  region  where  p  ^  has  gone  to 
zero.  The  point  P^  shown  in  the  figure  also 
lies  in  the  nonabsorbing  region  on  a  section 
of  the  ray  which  has  not  yet  entered  the 
absorbing  interior. 

The  intensity  at  a  point  P  is  then  given  by  the  integral 
T  (P. ) 

/V  '  • ' 

(p'»«p(T^  -T/)dT/  (2) 

T. 

If  has  been  determined  for  all  rays  going  through  a  given  point,  the  in¬ 
tegral  overall  directions 


/  ^  do  -  p  (B  -  I  )  d.o  (3) 

•/  dT  y  V  V 

is  the  difference  between  emitted  and  absorbed  power  per  unit  volume  and 
frequency  and  determines  the  rate  of  cooling  at  that  frequency. 

In  principle,  carrying  out  integrals  like  (2)  is  quite  straightforward  but 
because  of  the  strong  frequency  dependence  of  p  extensive  calculations  at 
many  frequencies  are  involved.  The  difficulty  which  arises  is  that,  even  for 
frequencies  which  lie  close  together,  the  major  contributions  to  the  integral 
for  I  may  come  from  widely  separated  points  on  the  ray.  The  variation  of 
the  point  of  origin  (i.  e,  where  the  integrand  in  equation  (2)  is  large)  with 
frequency  is  particularly  severe  for  transitions  between  bound  states  because 
they  lead  to  spectral  lines  where  the  absorption  coefficient  is  veryi^uch 
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larger  at  the  center  than  In  the  wings. 

Until  a  more  satisfactory  method  has  been  found,  one  is  forced  to  use 

approximate  procedures.  It  is  customary  to  use  smoothed  coefficients  for 

the  integration  of  the  transport  equation  in  which  the  contribution  due  to  lines 

is  distributed  in  frequency  as  if  they  were  broad  enough  to  overlap. 

The  current  status  of  calculations  to  obtain  the  absorption  coefficient  of 

hot  air  is  discussed  in  considerable  detail  in  a  forthcoming  article  by  Meyerott 
2 

and  Armstrong  .  Results  of  such  calculations  in  the  temperature  range  from 

1,  000°  K  to  12,  000°  K  and  for  density  ratios  relative  to  sea  level  extending 

from  10  to  10~  have  been  calculated  and  are  tabulated  as  a  function  of  wave- 

3 

length  between  1,  200  A  and  20,  000  A  .  The  actual  wavelengths  appearing  in 
the  tables  are  those  corresponding  to  values  of  hv  extending  from  0.  625  to 
10.  625  ev  in  steps  of  0.  25  ev. 

In  a  higher  temperature  range  extending  from  22,  000°K  to  220,  000°K  and 

in  nearly  the  same  density  range  as  before,  absorption  coefficients  have  been 

4 

reported  by  Armstrong,  Holland,  and  Meyerott  .  The  values  of  hv  extend 
from  1  ev  to  1,  000  ev,  i.  e.  ,  beyond  the  binding  energy  of  the  last  K  electron 
in  Owhich  is  871,  1  ev, 

A  number  of  assumptions  entering  into  the  calculations  of  absorption 
coefficients,  such  as  the  use  of  the  Born-Oppenheimer  approximation,  un¬ 
certain  f-values,  and  the  smoothing  out  prescriptions,  make  the  use  of  these 
coefficients  somewhat  suspect.  One  way  of  testing  them  for  the  conditions 
where  they  are  intended  to  be  used  is  the  calculation  of  the  spectral  distribu¬ 
tion  emerging  from  a  nuclear  fireball.  Such  a  calculation  has  been  carried 
out  for  a  1-KT  airburst  with  temperature  and  density  profiles  taken  from  the 
report  by  Brode  and  Meyerott  (ref,  RM-1851).  The  results  were  compared 
with  experimental  spectral  distributions  observed  from  such  small  airbursts. 
There  is  good  agreement  on  features  like  a  sharp  peak  at  5,  000  A  and  a 
sudden  ultraviolet  cutoff.  By  contrast,  the  computed  spectral  irradiances 
are  significantly  higher  than  the  experimental  values  in  the  infrared  region. 
This  would  result  if  the  free-free  absorption  coefficients  had  been  assumed 
too  low.  Recent  theoretical  and  experimental  results  confirm  that  the  free- 
free  coefficients  are  indeed  higher  than  the  ones  used  in  the  calculation. 
Between  the  UV  cutoff  and  the  maximum  of  the  observed  spectral  distribution, 
one  finds  a  number  of  discrete  spectra.  Prominent  among  these  are  the  O^ 
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Schumann  Runge  lines  which  appear  in  absorption.  Other  spectra,  such  as 

-  )  appear  originally  in  absorption  and  change  some  time  after  the  2nd 
maximum  to  emission.  The  Schumann  Runge  bands  are  of  particular  interest 
because  of  a  very  strong  sensitivity  to  temperature.  Figure  1  shows  the 
potential  energy  curves  of  the  two  electronic  states  of  0.  which  give  rise  to 
the  S-R  band  system.  Transitions  from  low-lying  levels  of  the  X  2)  g 
state  are  severely  restricted  by  the  Franck-Condon  principle.  Levels  with 

ft 

vibrational  quantum  numbers  v  like  3,  4,  or  5  are  much  more  likely  to  make 
optical  transitions,  but  these  levels  are  not  populated  until  the  temperature  is 
high  enough.  The  relative  intensities  of  the  various  lines  can  be  used  to  cal¬ 
culate  the  population  of  the  various  levels.  One  can  then  attempt  to  interpret 
these  populations  in  terms  of  temperatures.  At  early  times  in  the  history  of 
an  explosion,  that  attempt  gives  rise  to  contradictory  results,  which  is  an 
indication  that  these  states  are  not  present  in  thermodynamic  equilibrium 
concentrations.  At  times  around  the  first  minimum  the  situation  changes  and 
from  then  on  the  various  levels  have  a  temperature-determined  population. 

The  calculations  reported  in  reference  3  show  that  the  Schumann- Runge 
transitions  are  the  dominant  contributor  to  the  absorption  coefficient  at  low 
temperatures.  At  sea  level  density,  this  temperature  range  extends  to  about 
6,  000°  K  and  at  10*^  times  s-1  density  to  about  3,  000°K,  at  which  point  the 
NO  bands  take  over.  It  is  important  to  note  that  there  is  no  significant  back¬ 
ground  of  continuous  transitions  up  to  temperatures  like  5,  000°K.  This  fact, 
together  with  the  appearance  of  bands  in  the  observed  spectra  from  explosions, 
indicates  that  line  effects  may  invalidate  the  use  of  smoothed  out  absorption 
coefficients  in  calculating  the  radiative  transfer  of  energy. 

The  Schumann- Runge  and  other  band  systems  consist  of  exceedingly  large 
numbers  of  lines  which  are  characterized  by  their  wavelength,  width,  and 
some  measure  of  their  intensity  such  as  the  value  of  p  ^  at  the  center  of  the 
line  or  the  integral  ^  d  ^  .  The  wave  lengths  can  be  calculated  for  given 
vibrational  and  rotational  quantum  numbers  of  the  upper  and  lower  state, 

I  ti  s 

v'  J  and  V  J  ,  with  the  aid  of  constants  given  by  Herzberg  .  The  line 
widths  are  determined  by  pressure  broadening  and  are  assumed  to  be  pro¬ 
portional  to  the  pressure  with  a  constant  of  proportionality  which  has  been 
calculated  by  Margenau.  The  integrated  absorption  coefficient 
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Figure  1.  Levels  for  Schumann- Runge  transitions 
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2)  The  electronic  f  value,  which  for  0.  (S-R)  is  assumed  to  be 
f  a  0.  048  as  given  by  Treanor  and  Wurster  . 

3)  The  Franck-Condon  factors  g^n  ,  as  given  by  Fraser,  Jarmain,  and 
Nlcholls^. 


4)  The  line  strengths  S  n  ,  as  quoted  on  page  208  of  reference  5.  The 
integral  given  in  equation  (^)'^and  the  maximum  center  of  the  line 

are  directly  proportional  to  each  other.  The  factor  depends  on  the  line  shape 
and  varies  like  the  line  width  Av  .  Figures  2  and  3  show  the  envelopes  of 
various  bands  at  1,  000°  K  and  at  3,  000°  K. 


Figure  4  shows  the  fractional  transmission  for  several  thicknesses  x  of 
air  of  a  number  of  lines  in  the  (5-5)  band.  It  is  important  to  notice  how  these 
lines  become  wider  with  increasing  thickness.  For  very  thin  samples  the 
absorption  u  vx  would  be  proportional  to  the  thickness  but  this  increase  di¬ 
minishes  as  the  absorption  becomes  larger.  Thus,  the  centers  of  lines  where 
the  absorption  is  the  largest  grow  less  strongly  than  the  wings,  and  this  makes 
the  lines  become  wider.  Since  the  absorption  at  the  line  center  approaches 
its  full  value  in  a  relatively  thin  sample,  the  major  contribution  to  the  trans¬ 
mission  comes  from  the  wings  and  the  average  transmission  depends  critically 
on  the  minimum  of  ^ 

V 

The  procedure  of  integrating  the  transport  equation,  so  as  to  develop  as 
much  detail  as  shown  in  figure  4,  is  quite  laborious  and  with  more  than 
10,  000  and  possibly  as  many  as  100,  000  lines  the  amount  of  computation  would 
be  prohibitive.  There  are,  however,  certain  regularities  which  can  be  readily 
noticed  on  inspection  of  band  spectra,  which  suggest  a  simpler  modification 
of  the  above  procedure.  Figure  5  shows  the  densitometer  trace  of  a  typical 
band  spectrum  attained  with  a  high  resolution  Meinel  spectrograph.  The 
point  which  we  want  to  make  is  that  the  average  absorption  changes  quite 
slowly  with  wave  length.  The  fact  that  figure  5  is  an  emission  spectrum  does 
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Figure  2.  Intensity  distribution  in  Schumann-Runge  bands  at  T  =  1.  000°K 
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Figure  3.  Intensity  distribution 
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Figure  4.  Fr.ctlon.1  , 
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not  affect  this  argument  since  the  ratio  of  emission  and  absorption  is  also  a 
slowly  changing  function  of  wave  length.  The  gradual  change  of  absorption 
suggests  that  it  may  be  sufficient  to  solve  the  transport  equation  properly  in 
a  small  number  of  fairly  narrow  frequency  intervals  such  as  the  one  shown 
in  figure  4  and  use  these  results  to  obtain  averages  for  the  entire  spectrum. 

A  machine  code  to  do  this  has  been  written  and  has  been  very  nearly 
checked  out.  Details,  such  as  optimal  values  for  the  width  and  the  nuniber  of 
intervals,  have  to  be  found  by  experimenting  with  this  code.  It  is  estimated 
that  an  interval  may  contain  an  average  of  20  and  certainly  less  than  200  lines, 
and  that  there  may  be  a  total  of  20  such  intervals. 

The  general  procedure  then  is  as  follows.  For  a  given  time  in  the  history 


of  the  explosion  at  which  density  and  temperature  profiles  are  known,  one  can 
obtain  densities  and  temperatures  along  any  ray  through  the  fireball.  One  can 
therefore  calculate  absorption  coefficients  along  such  rays  for  the  selected 
frequencies  within  each  one  of  the  narrow  intervals.  By  integrating  the  trans¬ 


port  equation,  one  thus  obtains  the  intensities  I  which  are  then  averaged  for 

—  I  z'  ^ 

each  narrow  interval  to  form  I  ■ -  /  I  dv.  The  next  step  is  to  inter- 

Av 

polate  with  respect  to  frequency,  and  some  work  is  currently  underway  to 
devise  a  good  procedure  for  doing  this.  The  interpolated  I^is  then  integrated 
over  the  entire  spectrum  to  form 


00 


I  d  .  One  can  now  calculate  the  cooling  rate 

V  V  * 


this  into  the  energy  equation,  and  advance  the  hydrodynamic  calculation. 
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NUMERICAL  METHODS 

FOR  HYDRODYNAMICS  WITH  RADIATION  TRANSPORT 

Dr.  R.  A.  Grandey 

Aeronutronic  Division,  Ford  Motor  Company 

1.  THE  EQUATIONS  OF  RADIATION  TRANSPORT. 

The  integrodifferential  equations  of  radiation  transport  are  obtained  by 
simultaneously  solving  the  time  dependent  Boltzman  transport  equation  for 
the  radiation,  which  is  a  continuity  equation,  and  the  equation  expressing 
overall  conservation  of  material  plus  radiant  energy.  The  energy  equation 
is  written  in  the  form 

+  P  ^  ,  O. 

Here  E^  is  the  material  energy  density,  E^  the  radiation  density,  P  the  total 
pressure,  V  the  specific  volume  and  F  the  net  flux  of  radiation.  The  third 
and  fourth  terms  of  this  equation  are  expressed  as  integrals,  over  frequency 
and  direction,  of  the  Boltzman  equation. 

The  Boltzman  transport  equation  gives  the  time  rate  of  change  of  intensity 
of  radiation  of  frequency  v  traveling  in  a  direction  w.  This  intensity  will  be 
a  function  of  time  and  position; 


where  r  is  the  position  vector  and  t  the  time.  The  transport  equation  may 
be  written  generally  in  the  following  form: 

+  (w-  7)  =  -  k*’  e  *  pv 

The  left  side  of  this  equation  represents  1/c  times  the  total  time  derivation  of 

I  The  first  term  on  the  right  represents  photons  removed  by  absorption 

^  ^  (a) 

processes,  being  the  absorption  coefficient.  The  second  term  on  the 

right  represents  photon  production  by  emission  from  the  material,  the  third 
term  those  photons  lost  by  scattering,  and  the  fourth  those  added  by  scattering 
from  other  directions,  "W  ,  into  the  direction "w.  In  this  equation  P  is  the 


cli 
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is) 

material  density,  ^ the  emission  coefficient  and  ,  scattering  coefficient. 


(2) 


Our  problems  will  assume  slab  or  spherical  symmetry,  in  which  case 
I  is  a  function  of  only  one  coordinate  x  and  the  angle  0  between  the  di- 
rection  of  symmetry  and  the  direction  of  travel  of  the  photon  in  question: 

^y,W  ""  j 

where  p,  s  cos  0.  We  may  write  equation  (1)  as 

+  j P  lyy 

K  is  the  total  "absorption"  coefficient,  i.  e.  ,  the  sum  of  the  pure  ab- 

^  I 

sorption  and  pure  scattering  coefficients.  oi"  =  0  or  2  for  slab  or  spherical 
geometry  respectively. 

The  net  monochromatic  flux,  ,  traversing  unit  area  is  a  vector  in  the 
x-direction  with  magnitude  2xc  ^jidp.  The  total  flux  F  is  defined  by 

the  y/* F.  dv  or:  "• 


F-  F,  -  a»c  j 


(3) 


O  -I 

The  radiation  energy  density  is  the  integral  of  ^  over  direction  and 
frequency,  i.  e. 

£’r*2»/7lv,^4,d»» 

o  -I 

Thus  the  energy  equation  may  be  written  as: 

w  A 


(4) 


(5) 


140 


TN-61-29 


The  transport  equation  is  applied  to  assist  in  evaluating  the  integrals  appear¬ 
ing  in  this  equation. 

The  term  ^  transport  equation  represents  both  spontaneous  and 

induced  emission.  Its  form  is  sensitive  to  whether  or  not  local  thermodynamic 
equilibrium  exists  for  the  material.  If  local  equilibrium  is  assumed,  it  may 
be  expressed  in  terms  of  the  absorption  coefficient  and  the  Planck 

spectrum.  By  combining  this  with  the  absorption  term,  equation  (2)  becomes 

1  4,  u 

t  at  ^  ax  ax 


In  our  problem,  scattering  may  often  be  ignored  in  comparison  with  the 
emission  process.  In  this  case. 


Generally,  we  must  solve  equations  (5)  and  (7)  subject  to  the  appropriate 
boundary  conditions.  Introduce  an  integrated  intensity 


CO 
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by  integrating  equation  (7)  over  dv: 


ti  ®® 

aT  _  _  f  \-'t 


(8) 


dv 


where  X^,  the  transport  mean  free  path,  ie  defined  by  ^-1  dy  dv 


and  a  is  the  Stefan-Boltzman  constant.  The  second  term  in  equation  (S)  may 
be  eliminated  by  integrating  (8)  over  d[i: 

i  it  {l  dv  -  T  Jj'  hj  4«  dv  - 

K'^T dji  dv  -|x  I  ^ 

Our  problem  Is  to  solve  equations  (8)  and  (9). 

All  equations  are  linear  in  the  radiation  intensities.  Let  us  therefore  write 

the  monochromatic  intensity  as  the  sum  of  two  terms:  one  including  radiation 

due  to  material  at  a  local  temperature  T.  and  the  other  radiation  coming  from 

an  external  source  in  transit  through  the  material.  The  first  term  will  include 

radiation  flowing  into  the  material  from  neighboring  material  not  necessarily 

p  •  E 

at  exactly  the  same  temperature.  That  is,  define  I  si  -t-  I  ; 

I  ■  I  -t-  I  ,  where  the  superscripts  stand  for  Planck  and  external,  and  write 
separate  continuity  equations  for  the  two  components: 


(10) 
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£ 

KA 


(11) 


Specializing  to  the  case  of  ^  s  0  for  simplicity,  the  energy  equation  (9)  may 
be  written  as: 


A;'dT*-  jt  V;  i;,  4.  jvI 


(12) 


-fl' 

-*0  '~\ 


X  I 


c/y  =  o . 


The  diffusion  approximation  for  the  Planck  spectrum  is  obtained  by  assuming 
that  pK  is  very  large,  or  that  \  is  small  compared  to  dimensions  over 

*  p 

which  quantities  vary  appreciably.  In  this  case,  an  expansion  of  I  as  a 

power  series  in  |i  converges  rapidly.  The  result  of  such  an  approximation  is 

r* 

the  usual  diffusion  equation;  plus  an  additional  term  involving  I  : 

oo  \  Q3) 

-|  Ix^Ar  -c jy  J^Jy  =  O 

-I 


where  X^,  the  Rosseland  mean  free  path,  is  defined  by 

’  It 

In  several  applications,  it  may  be  possible  to  make  use  of  equation  (13)  in 
treating  the  Planck  radiation,  using  equation  (11)  to  treat  the  transport  of 
external  radiation  through  the  material.  In  cases  where  the  diffusion  approxi¬ 
mation  is  not  valid  for  the  Planck  spectrum,  a  second  appreciable  simplifi¬ 
cation  is  still  possible.  Returning  to  equation  (7),  the  right  hand  side  of  this 
equation  may  be  written  X  (J  -I  -I  ).  The  integral  of  the  first  term 

V  V  V I  V  I  [1 
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of  this  expression  led  to  introduction  of  the  transport  mean  free  path,  .  It 

is  seen  to  be  appropriate  to  use  this  same  mean  free  path  as  an  approximation 

P 

in  equations  (10)  and  (12)  in  the  term  involving  I  : 


V. 


in  time  as  it 


We  shall  have  to  follow  the  changing  frequency  spectrum  of  i 
transports  through  the  material.  The  number  of  frequency  groups  necessary 
and  the  expense  of  treatment  will  depend  on  the  shape  of  the  spectrum  and  the 
precise  dependence  of  X  ^  on  v . 

Inclusion  of  the  pure  scattering  terms  and  the  extra  terms  for  the  spher> 
ically  symmetric  case  do  not  introduce  any  difficulties  and  have  been  omitted 
only  for  convenience  in  the  discussion. 


2.  NUMERICAL  INTEGRATION  OF  THE  TRANSPORT  EQUATIONS. 

Three  approaches  have  been  taken  toward  solving  the  transport  equations. 
The  chief  problem  is  treating  the  angular  integration.  An  immediate  extension 
of  the  diffusion  approximation  is  to  use  higher  order  terms  in  p,  or,  equiva¬ 
lently,  to  expand  I  (p)  in  Legendre  polynomials,  since  one  can  then  make  use 
of  the  various  recursion  and  orthogonality  relations  existing.  This  approach 
leads  to  a  set  of  L  1  equations  for  the  Legendre  coefficients  if  polynomials 
up  to  order  L  are  used  in  the  expansion.  I  may  be  defined  by  an  L  -**  1  com¬ 
ponent  vector  consisting  of  the  Legendre  coefficients  and  finite  difference 
methods  applied  to  the  resultant  vector  equation.  The  coefficients  in  the  vector 
equation  are  L  1  by  L  1  matrices.  The  difficulty  of  this  so-called 
"Spherical  Harmonic  Method"  is  that  an  implicit  difference  scheme  is  required 
for  unconditional  stability  and  the  scheme  must  be  implicit  both  with  respect 
to  the  space  coordinate  and  the  "1  coordinate,  "  yielding  very  cumbersome 
equations. 

A  second  approach  is  based  on  the  "method  of  discrete  ordinates."  The 
distribution  function  is  sampled  in  each  of  a  finite  number  of  directions,  and 
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separate  equations  written  for  the  individual  directions.  The  two  approaches 
are  related,  inasmuch  as  if  L  +  1  discrete  angles  are  used  for  sampling,  the 
solution  can  be  fitted  with  the  first  L  spherical  harmonics.  Wick  and 
Chandrasekhar  developed  this  approach  independently  using  a  Gaussian  quadra¬ 
ture  approximation  for  the  angle  integrals.  It  can  be  shown  that  the  results 
from  this  approach  are  identical  to  results  of  the  spherical  harmonic  method. 

This  method,  which  uncouples  the  equations  angularly,  is  much  more 
tractable  than  the  spherical  harmonic  method  if  it  is  desired  to  obtain  un¬ 
conditionally  stable  finite  difference  equations. 

3.  THE  CARLSON  S  METHOD. 

n 

The  method  is  also  a  discrete  ordinate  scheme  in  carrying  out  the 
integration  over  angle.  It  differs  from  the  Wick-Chandrasekhar  (W-C)  method 
in  several  respects.  The  centering  with  respect  to  p  is  done  differently  and  a 
trapezoidal  rather  than  Gaussian  quadrature  is  used  to  perform  the  integration 
over  angle. 

The  S  method  divides  the  p  interval  into  n  subintervals;  -1  »  <p,  <-  - 

n  2^  ° 

--<u  a+1.  Ifp,-  -1+  /  ,  the  p-mesh  is  said  to  be  standard.  Other 

^  n  i  '  n 

choices  may  be  appropriate  to  some  problems. 

In  the  W-C  method,  the  transport  equations  are  centered  at  the  p-mesh 
points.  S  centers  these  equations  at  the  midpoint  of  the  p  -intervals: 

X  i*'/r.  —  1/  ^  r  / » 7x  U  ~^(i*  ’/%  ]  ^  r  I  (16) 


This  set  of  n  equations  in  n  +  1  unknowns  is  completed  by  adding  the 
equation  centered  at  1  z  o  (pz  -1): 
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Assume  that  i  may  be  represented  by  K  frequency  groups 
B  K 

T  ^  y  1 

xy  u  Lj  write  separate  transport  equations  analagous  to 


equations  (16)  and  (16a)  for  the  K  components: 


C  iX  2X  ^  IjV/s 


(17) 


t  3+  dx  l»/» 


(17a) 


The  number  of  angular  groups  required  for  accuracy  may  be  different  for 
P  E 

I  and  I  .  It  is  planned  to  allow  for  using  an  integration  for  the  Planck 
spectrum  and  an  Sj  integration  for  the  "external"  radiation. 


3.  1  The  Difference  Equations. 

Equations  (17)  and  (17a)  will  be  used  to  illustrate  Carlson's  differencing 
methods.  A  more  complicated  procedure  is  required  if  anisotropic  scattering 
(like  Thomson  scattering)  is  included  but  we  shall  not  discuss  this  in  the 
interests  of  simplicity  and  time.  We  shall  drop  the  K  subscript  for  simplicity. 

Equation  (17)  is  written  in  terms  of  the  intensities  at  p^mesh  points  j  and 
j+  Iby  performing  an  integration  over  p  under  the  assumption  that  the  intensity 
I  is  linearly  dependent  on  p  in  the  interval  (a  trapezoidal  approximation): 


^Llf  s  -f 


Aj  ~Aj.t 


aI/,  T  -  T- 

AJ  -Ah 
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«  Wt  it  '  -Mi; 

The  above  are  introduced  into  equation  (17)  and  an  integration  performed  over 
^  from  to  Hj.  We  note  that 

J-l 


where: 


-  i  Uj4i  t/Tj-, ) 

~  3  *  ^Mj; ) 

b  j  s  z  (3  Mi  -MiMj.,  vYa#  ) 

^  {Mi -Mi;) 

Introducing  these  into  equation  (17)  and  multiplying  by  2/L: 


(18) 


(17a)  simply  becomes: 


[c  It  lx  'Jlo 


(18a) 
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Normally,  in  spherical  systems  the  radi'ition  streams  in  the  direction  of 
increasing  |j. .  The  angular  dependence  is  treated  by  arguing  that  the  differ¬ 
encing  should  be  backward  with  respect  to  the  direction  of  flow.  It  is  on  this 
consideration  that  equations  (16a)  and  (17a)  are  for  p  5  -1  rather  than  for 
p  ■  +  1.  Equation  (18a)  is  solved  for  which  is  inserted  into  the  equation 

involving  and  which  is  solved  for  etc. 

We  shall  now  introduce  spatial  coordinates  which  define  a  spatial  mesh 
with  an  interval  A  =  X  -X  ,  and  the  temporal  coordinates  t'^which  define  a 
temporal  mesh  with  an  interval  =  t‘  -t'^  .  Let  us  examine  equations  (18) 
and  (18a)  centered  at  time  t  and  at  the  spatial  coordinate  X 

m 


a  distance  ac&  from  point  e.  They  traveled  along  path  1,  2,  3,  or  4  depending 

on  the  value  of  ac&.  The  spatial  and  temporal  derivatives  occurring  in  equation 

(18)  are  formed  along  the  appropriate  legs  of  the  triangle  traversed  by  the 

photons  in  question.  For  example,  if  the  photons  follow  path  1  (if  acb>A  ), 
d  I 

the  j  appearing  inJ18)  is  differenced  as 

*  Ayyy 
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Generally,  similar  equations  will  exist  for  each 


The  result  of  the  above  will  be  an  implicit  finite  difference  equation  for 

I,  ^  X  )  =  Generally,  similar  equations  will  exist  for  each 

J  A 

frequency  group,  and  it  will  be  necessary  to  define  1  <7  j 

jt 

It  should  be  noted  that  the  transport  equations  are  coupled  angularly  through 
the  scattering  integral  and,  in  the  case  of  spherical  symmetry,  through  the 
term  ^ .  In  slab  geometry  with  no  scattering  the  equations  are  not 
coupled  angularly  and  an  -type  calculation  would  not  be  required, 

A  third  approach  was  mentioned.  It  is  also  possible  to  write  the  angular 


term 


dependence  in  an  integral  form  and  apply  iterative  techniques  which  will  con¬ 
verge  to  the  correct  solution^. 
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SHOCK-INDUCED  REACTIONS* 

Dr.  O.  A,  Nance 
E.  H.  Plesset  Associates,  Inc. 

In  discussing  the  possibility  of  shock-induced  reactions  in  condensed 
phases,  we  find  that  even  reliance  on  the  weak  reed  of  chemical  reaction 
theory  eludes  us.  We  can  only  hope  to  avoid  the  nonsequitur  that  experi¬ 
mental  evidence  provides  many  examples  of  this  phenomenon. 

The  framework  of  gas  phase  reaction  kinetics  is  dually  supported,  by 
thermodynamics  and  by  the  "absolute"  rate  theory.  Classically,  reactions 
are  supposed  to  occur  by  "activation"  of  a  molecule  through  collision  or  ab¬ 
sorption  of  radiation.  If  the  activated  molecule  simultaneously  or  subse¬ 
quently  encounters  another  appropriate  molecule,  a  bimolecular  reaction 
may  occur.  Alternatively,  the  activated  molecule  may  decompose  in  a  uni- 
molecular  reaction.  To  various  degrees  of  sophistication,  we  may  calculate 
the  number  of  collisions  occurring  as  a  function  of  concentration  and  tempera¬ 
ture. 

Calculation  of  the  rate  of  reaction  from  the  number  of  collisions  alone 
usually  leads  to  a  gross  overestimate  of  the  rate.  The  next  level  of  approxi¬ 
mation  requires  an  estimate  of  the  energy  involved  in  collisions  and  an  esti¬ 
mate  of  the  probability  of  transferring  an  effective  amount  of  energy  to  one 
of  the  colliding  particles.  The  first  formal  expression  is  due  to  Arrhenius 
and  takes  the  form  k  =  Ae"  *  where  k  is  the  specific  rate  constant,  A 

is  a  frequency  factor  and  E  is  the  energy  of  activation.  This  expression 
was  inherently  an  empirical  description  of  the  variation  of  reaction  rate  with 
temperature.  "A"  can  be  related  to  the  number  of  collisions  but  with  only 
indifferent  success  except  in  the  case  of  simple  bimolecular  reactions  be¬ 
tween  atoms. 

If  we  add  the  concept  of  the  Maxwell  distribution,  we  gain  the  concept  of 
the  fraction  of  molecules  capable  of  undergoing  effective  collisions  and  the 
concept  of  at  least  a  quasi-equilibrium  state  which  will  maintain  the  popula¬ 
tion  of  "hot"  molecules.  Use  of  this  new  information  improved  the  calculation 

*  This  report  is  based  on  preliminary  results  of  work  sponsored  by  the 
Richfield  Oil  Corporation.  Potential  applications  to  petroleum  technology 
are  subject  to  patent  rights  sought  by  that  Company. 
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of  the  frequency  factor  in  simple  cases  but  gave  no  clue  to  the  predictability 
of  the  energy  of  activation. 

With  the  advent  of  the  absolute  rate  theory,  a  major  practical  advance 
occurred;  unfortunately  the  quantitative  gains  and  theoretical  validity  were 
severely  limited.  In  this  theory  we  consider  a  potential  energy  surface  rep¬ 
resenting  the  reacting  system  in  the  pertinent  region  of  configuration  and 
momentum  space.  One  assumes  a  quasi-equilibrium  distribution  of  the  macro 
state  with  respect  to  this  surface  and  examines  those  portions  of  the  potential 
energy  surface  with  which  externally  recognizable  chemical  changes  maybe 
correlated.  In  the  customary  notation,  we  write 

^AS  -AH  VrT  ^  kT  Q^^-AE8/RT 


Where  is  the  specific  rate  constant 
T  is  absolue  temperature 
h  is  Planck's  conscant 

AE§  is  the  difference  in  "ground  state"  between  normal  and  activated 
species 

Q  represents  a  partition  function 

AS^iS  the  entropy  of  activation 

AH^  represents  part  of  the  energy  of  activation 


The  superscript  ^  refers  to  the  activated  species  and  the  subscript  f 
refers  to  the  quantity  after  the  "reaction  coordinate"  is  removed.  The  pro¬ 
duct  of  Qj  over  "i"  is  necessary  when  several  reactant  species  are  involved. 
It  will  be  necessary  to  refer  the  reader  to  standard  texts  or  reference  works 
for  discussion  of  the  development  of  these  expressions. 

The  point  of  outlining  the  background  briefly  here  is  to  permit  a  dis¬ 
cussion  of  the  situations  in  which  misapplication  is  possible.  The  first  set 
of  terms  defines  the  rate  in  thermodynamic  variables;  this  is  of  value  pri¬ 
marily  in  recognition  of  the  source  of  a  temperature  dependent  and  a  tem¬ 
perature  independent  term,  together  with  clues  to  the  application  of 

kT 

statistical  mechanical  concepts.  One  will  note  that  the  factor  appears 
in  both  expressions  and  that  it  has  been  borrowed  from  AH  ^  in  one  case 
and  from  in  the  other  (noted  by  demotion  of  the  symbol  /  to  a  subscript). 
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In  the  second  expression,  the  partition  functions  are  taken  with  respect  to 

a  ground  state  in  each  case  and  the  ground  state  factors  are  combined  in  the 

term  e  In  addition,  the  partition  function  has  been  depleted  by  the 

factor  kT;  this  factor  is  usually  interpreted  as  the  partition  function  of  the 
■E~ 

"reaction  coordinate",  treated  as  a  very  "stiff"  vibration.  That  is,  it  is 
assumed  that  an  important  characteristic  of  the  activated  complex  is  the 
occurrence  of  an  appropriate  amount  of  energy  in  a  "degree  of  freedom"  of 
the  complex  such  that  motion  in  that  one  coordinate  can  produce  the  observed 
macroscopic  reaction.  In  terms  of  the  potential  energy  surface,  motion  in 
this  coordinate  leads  from  "reactant"  to  "product"  portions  of  the  energy  sur¬ 
face.  Roughly  speaking,  AEq  ,  is  the  least  upper  bound  of  the  energy  required 
to  surmount  the  potential  barrier  and  thus  is  in  excess  of  the  quantities  AH^ 
or  E  which  are  averages  of  the  energy  barrier  from  a  given  level  to  the 
transition  level,  weighted  by  the  energy  distribution  for  the  reactant  species. 

For  gas  phase  reactions  of  simple  species,  the  partition  functions  may  be 
approximated  and  estimates  may  be  made  of  the  partition  functions  for  postu¬ 
lated  activation  configurations.  Furthermore,  we  have  a  reasonable  idea  of 
the  use  of  temperature,  both  to  describe  the  energy  partition  and  to  calculate 
the  number  and  effect  of  collisions.  We  can  also  estimate  the  effect  of  pressure 
on  the  rate  of  reaction  and  on  the  distribution  of  products,  if  several  reaction 
paths  are  possible.  Increasing  the  pressure  on  a  gas  raises  the  temperature 
and  increases  the  concentration  of  reactant  species;  it  also  increases  the 
probability  of  deactivation  and.  unless  the  activated  complex  is  more  compact 
than  the  separated  reactants,  may  lead  to  increase  or  decrease  in  reaction 
rate,  depending  on  the  specific  reaction  mechanism. 

In  the  gas  phase,  the  probable  collision  c*rientations  and  transfer  of  energy 
between  degrees  of  freedom  is  also  somewhat  amenable  to  treatment.  In 
particular,  the  energy  content  is  not  readily  exchanged  between  translation 
and  either  rotation  or  vibration,  i.  e.  .  there  is  difficulty  in  tran.sfer ring  * 
translational  energy  and  there  is  a  corresponding  persistence  of  vibrational 
excitation  when  it  does  occur.  The  effectiveness  of  vibrational  energy  in 
complex  molecules  is  complicated  as  well  by  the  possibility  of  delocalization 
and  localization  which  arises  from  enharmonic  contributions  to  the  potential 
functions. 
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Statlitical  theory  of  tranaport, 

'  The  fundamental  questions  which  apply  here  simply  emphasize  the  limita¬ 

tions  of  reaction  theory  in  general.  One  of  the  fundamental  difficulties  here- 
/  tofore  has  been  that  theories  predicting  reaction  rates  have  postulated  equili¬ 

brium  states  and  used  thermodynamic  expressions  more  from  desperation 
than  sagacity.  A  fruitful  and  logical  point  of  attack  must  lie  in  the  larger 
question  of  description  of  the  approach  to  equilibrium  and  irreversible  pro¬ 
cesses  in  statistical  mechanical  terms. 

In  attempting  to  predict  kinetic  processes  we  are  attempting  to  assess 
the  behavior  of  perturbed  regions,  molecular  aggregates,  or  individual 
molecules  in  a  macroscopic  system  which  may  be  at  or  near  equilibrium. 

This  point  has  long  been  recognized,  e.  g.  ,  Kirkwood  made  significant  contri¬ 
butions  and  was,  at  the  time  of  his  death,  a  leader  in  this  technique  of 
approach. 

It  was  the  author's  original  intent  to  discuss  nonequilibrium  theory  in 

some  detail  but  it  became  increasingly  obvious  that  such  a  discussion  was 

impractical  under  the  circumstances.  The  reader  who  wishes  to  pursue  the 

matter  is  referred  to  an  excellent  review  article  by  S.  A.  Rice  and  H.  L. 

Frisch  which  appeared  in  vol.  11  of  the  Annual  Review  of  Physical  Chemistry 

(pps.  187  -  272)  (Ann.  Reviews,  Inc.,  Palo  Alto,  Calif.  -  1960)  or  to  recent 

/ 

papers  by  Kirkwood  and  his  collaborators,  particularly  S.  A.  Rice. 

Rcftctiong  ifl  gondgna^d-phaasg. 

With  the  foregoing  discussion  in  mind,  great  audacity  is  required  in  ap¬ 
plying  the  theory  of  gas  phase  reactions  to  condensed  phases  and  shock  phen¬ 
omena.  To  the  extent  that  shocks  can  raise  the  local  energy  content  sub¬ 
stantially,  we  can  begin  to  consider  reactions  in  terms  of  the  general  theory. 
Activation  leading  to  unimolecular  reactions  can  probably  be  accomplished  by 
collisions  between  neighboring  atoms  or  molecular  groups.  In  the  sonic 
approximation,  the  transfer  of  energy  to  vibrational  modes  of  a  crystalline 
structure  is  known  to'be  an  effective  means  of  energy  transfer  and  we  would 
expect  individual  species  to  be  subjected  to  periodic  disturbances.  Since  these 
disturbances  have  periodscharacteristic  of  the  lattice  modes,  the  coupling  to 
^  atomic  vibrations  of  much  higher  frequency  would  be  of  comparatively  low 

efficiency.  In  the  transonic  and  strong  shock  regimes,  the  efficiency  of 
transfer  should  appreciably  increase.  Furthermore,  loss  of  bond  vibrational 
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excitation  should  be  inhibited  by  the  factors  which  by- passed  molecular 

excitation  in  the  sonic  case.  Some  measure  of  the  vibrational  excitation  can 

be  obtained  in  gases  by  observing  the  dispersion  of  supersonic  waves;  there 

it  is  observed  that  vibrational  excitation  by  one  or  more  quanta  can  be  ob- 

2  5 

tained  in  the  range  of  10  -  10  sound-induced  collisions.  In  solids,  the  situ¬ 
ation  is  more  complicated,  theoretically  and  experimentally,  but  it  has  been 

_9 

observed  that  the  lifetime  of  vibrational  excitation  is  of  the  order  of  10 
seconds  in  many  cases. 

The  partition  functions  for  an  organic  material  should  be  markedly  differ¬ 
ent  in  gaseous  and  condensed  phases  because  in  condensed  phases: 

1.  the  translational  modes  are  strongly  coupled  to  the  lattice  modes; 

2.  the  rotational  degrees  of  freedom  are  severely  restricted  by  the 
limited  free  space  in  condensed  phases; 

3.  the  vibrational  modes  of  molecules  probably  have  higher  effective 
anharmonicity  because  of  the  intermolecular  force  fields,  but  are 
sufficiently  different  in  frequency  from  the  lattice  modes  to  have  appre¬ 
ciably  different  energy  content; 

4.  there  is  a  high  degree  of  degeneracy  in  the  lattice  modes. 

The  entropy  of  activation  in  the  gas  phase  may  be  considered  in  two  parts. 
The  first  of  these  arises  from  the  relative  orientation  of  the  colliding  mole¬ 
cules  at  and  during  collision,  and  the  second  derives  from  configurations 
assumed  during  the  free  oscillation  and  rotation  subsequent  to  the  collision. 

In  the  case  of  condensed  phases,  both  sources  of  variation  of  orientation  are 
severely  restricted  by  the  nature  of  the  state,  particularly  in  cases  where 
mechanical  treatment  or  preferred  bonding  configurations  are  present  during 
reaction  or  fabrication  of  the  material. 

By  analogy  with  phenomena  observed  following  irradiation  of  molecular 
aggregates,  the  predominant  reactions  appear  to  be  unimolecular  dissociation 
and  reaction  of  activated  sites  with  neighboring  atoms  or  groups.  A  special 
case  is  that  in  which  hydrogen  atoms,  produced  by  bond  rupture,  combine 
locally  or  diffuse  to  react  at  sites  somewhat  distant  (usually  100  angstroms 
or  less)  from  the  point  of  origin. 

CgnYcntioaal  rcftctiong- 

With  the  differences  noted,  it  should  be  possible  to  treat  reactions  in 
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condensed  phases  by  the  same  broad  theoretical  treatment  and  some  tech¬ 
niques  parallel  to  those  used  in  gas  phase  kinetics.  The  temperature  is  now 
replaced  by  the  internal  energy  of  the  material  (with  due  regard  to  potential 
difference  in  partition  functions);  the  expected  duration  of  the  reactions  is 
that  of  the  pulse  duration  plus  one  or  more  characteristic  relaxation  times 
which  might  be  distinguishable  by  variation  in  product  distribution. 

The  orientation  factors  are  fairly  stable  and  might  be  treated  either  by 
straightforward  probability  theory  for  isotropic  material,  or  by  considerations 
of  random  orientation  plus  the  concept  of  crystallites  familiar  to  workers  in 
the  field  of  plastics.  Since  some  control  over  samples  is  possible,  one  might 
be  able  to  use  variation  in  sample  fabrication  to  vary  reaction  rate  and  mech¬ 
anism  in  experimental  studies. 

Experimental  clues. 

The  facilities  required  for  study  of  shock-induced  reactions  in  condensed 
phases  are  not  readily  accessible.  One  requires  instrumental  bunkers  in 
locations  fairly  remote  from  living  areas,  a  supporting  staff  capable  of 
sophisticated  manipulation  of  high  explosives,  and  rather  elaborate  sample 
containers  which  can  survive  the  shocks  and  preserve  the  samples  for  post 
reaction  analysis.  Needless  to  say,  the  experiments  are  also  rather  expen¬ 
sive. 

Some  work  has  been  done  in  support  of  the  "Plowshare"  program  by  LRL 
and  some  work  has  been  supported  by  the  Richfield  Oil  Corporation  at  SRI. 
There  are  also  scattered  bits  of  data  which  come  as  by-products  of  military 
programs. 

One  can  summarize  the  available  data  crudely  as  follows: 

1)  severe  decomposition  of  organic  materials  is  produced  in  times  of  the 
order  of  microseconds  by  shocks  of  Z50  kb  or  higher; 

2)  limited  changes  in  materials  occur  in  the  same  time  intervals  down  to 
the  order  of  20  kb. 

Calculation  program. 

Preparation  of  a  numerical  machine  code  is  underway  in  support  of  this 
program.  The  code  consists  of  a  conventional  Lagrangian  hydrodynamics 
calculation  supplemented  by  auxiliary  sections  for  reaction  calculation.  For 
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materials  in  which  reaction  is  possible,  the  logic  of  the  hydrodynamic  calcu¬ 
lation  is  interrupted  at  the  point  where  the  pressure  and  internal  energy  have 
been  calculated. 

Two  possible  types  of  reaction  have  been  provided.  The  first  type  is  a 
''slow”  reaction  for  which  the  rate  equation  has  a  unimolecular  form  and  the 
rate  constant  is  written  in  superficial  Arrhenius  form,  i.  e.  ,  the  rate  is  first 
order  in  concentration  and  the  rate  constant  consists  of  the  product  of  a  con¬ 
stant  and  an  exponential  function  of  internal  energy.  During  each  time  cycle, 
the  differential  change  in  concentration  and  energy  is  calculated  for  the  pre¬ 
dicted  reaction  and  the  results  used  to  correct  the  hydrodynamic  quantities. 

In  the  event  that  changes  produced  by  the  reaction  are  comparable  to  changes 
produced  by  the  mechanical  factors  of  the  shock,  an  additional  iteration  is 
required.  In  the  cases  investigated  to  date,  the  energy  absorbed  or  produced 
by  the  reaction  is  a  small  percentage  of  the  total  and  no  additional  references 
to  the  hydrodynamic  portion  are  required.  It  is  obvious,  however,  that  highly 
endothermic  or  exothermic  reactions  could  materially  affect  shock  propaga¬ 
tion.  A  limiting  case  is  that  encountered  in  detonation  of  explosives,  in  which 
the  reaction  energy  is  the  mechanism  for  generation  of  the  shock.  It  may  be 
noted  in  passing  that  some  of  the  difficulties  in  the  theory  of  detonation  which 
have  been  ascribed  to  state  behavior  are  probably  due  to  neglect  of  the  effect 
of  high  pressure  on  the  reaction  rate. 

Some  preliminary  calculations  have  been  made  on  the  effect  of  an  under¬ 
ground  nuclear  explosion  on  petroleum  bearing  formations.  An  underground 
nuclear  explosion  is  characterized  by  initially  high  attenuation  of  peak  pressure 
with  distance,  because  of  both  geometric  and  energy  absorption  factors.  At 
the  same  time  it  is  characterized  by  persistence  of  pressure  which  is  greater 
by  several  orders  of  magnitude  than  the  impulse  duration  of  conventional  ex¬ 
periments  with  high  explosives.  Using  reaction  constants  derived  from  con¬ 
ventional  petroleum  technology,  the  calculations  indicate  extensive  chemical 
reaction  in  regions  near  the  detonation  with  a  sharp  decline  in  the  amount  of 
reaction  at  distances  where  the  peak  pressure  is  below  about  100  kilobars. 

This  result  is  to  be  expected  because  of  the  conservative  approach  in  assign¬ 
ing  reaction  parameters,  but  the  results  compare  well  with  the  limited  ex¬ 
perimental  data. 
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We  have  also  proposed  that  there  is  probably  another  type  of  reaction 
accompanying  shocks  and  which  we  have  designated  the  "fast"  reaction.  It 
has  been  observed  in  shock  tube  experiments  and  in  detonation  of  explosives 
that  the  observed  rate  of  reaction  is  frequently  higher  than  that  predicted  by 
the  theoretical  "temperature".  J.  von  Neuman  commented  on  this  point  in 
discussing  detonation  theory  of  explosives  and  remarked  that  the  behavior 
seemed  more  like  that  resulting  from  a  mechanical  blow  than  from  thermal 
excitation.  He  had  calculated  a  reaction  rate  by  a  method  comparable  to  our 
"slow"  reaction  and  found  the  predicted  rate  to  be  less  than  the  observed  rate 
by  a  factor  beyond  that  which  could  be  explained  by  parametric  ranges  avail¬ 
able  in  his  theory.  Unfortunately,  he  did  not  pursue  the  point. 

We  have  made  provision  for  this  second  type  of  reaction  in  the  code  but 
have  not  used  that  section,  pending  further  theoretical  work.  One  can,  how¬ 
ever,  discuss  some  possibilities  in  phenomenological  terms.  One  postulates 
that,  as  the  shock  strength  increases,  there  is  a  point  at  which  the  particle 
velocity  exceeds  that  for  which  the  predominant  energy  transfer  is  to  bulk 
modes  of  the  material  (lattice  modes  for  crystals  or  crystallites).  The  effect 
would  then  be  to  transfer  increasing  fractions  of  the  dissipated  shock  energy 
to  vibrations  of  the  molecules.  On  the  basis  of  dispersion  of  supersonic  waves, 
there  is  the  indication  that  bond  vibrational  excitation  persists  with  character- 

_9 

istic  half-life  of  the  order  of  10  second.  With  vibration  periods  of  the  order 
-13 

of  10  second,  there  is  am  pie  time  for  reaction  or  internal  transfer  of  the 
vibrational  excitation.  Since  unimolecular  decompositions  are  the  result  of 
localization  of  energy  in  vibrational  modes,  the  effect  of  relatively  direct 
transfer  of  shock  energy  to  vibrational  excitation  would  produce  reaction 
rates  far  higher  than  those  predicted  on  the  basis  of  the  bulk  "temperature" 
of  the  material.  At  present,  quantitative  prediction  of  the  effect  seems 
impossible,  and  consideration  of  secondary  factors  such  as  the  effect  of 
orientation  seems  pointless  except  for  speculation.  Nevertheless,  there  is 
a  possibility  that  recognition  of  the  qualitative  features  of  such  a  reaction 
might  permit  the  design  of  fruitful  experiments. 

From  the  standpoint  of  the  use  of  a  hydrodynamic  code,  these  considerations 
emphasize  the  limitations  of  a  formulation  which  uses  the  Richtmyer  -  von 
Neuman  "Q",  since  the  region  of  the  shock  front  which  is  so  grossly  dis¬ 
torted  by  the  use  of  an  artificial  viscosity  is  Just  the  region  in  which 
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accelerated  reactions  would  occur.  In  addition  to  consideration  of  experiments 
to  establish  reaction  parameters,  we  are  considering  known  and  pursuing  novel 
techniques  for  treating  the  shock  front  phenomena  in  condensed  phases.  At 
present,  we  have  no  significant  contribution  to  report. 
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CALCULATION  OF  SPALL  BASED  ON  A  ONE- DIMENSIONAL  MODEL 

Mark  Wilkins 

Lawrence  Radiation  Laboratory 

INTRODUCTION. 

Consider  a  one-dimensional  model  of  an  HE  metal  plate  system.  Since  the 
pressure  pulse  from  the  HE  at  the  HE- plate  interface  decays  rapidly  with  time, 
a  rarefaction  wave  proceeds  forward  into  the  plate,  degrading  the  advancing 
pulse  from  the  rear.  Upon  reaching  the  front  surface  of  the  plate  where  the 
normal  stress  must  be  zero,  the  pressure  pulse  is  reflected  as  a  receding 
rarefaction  wave.  This  wave  drops  the  pressure  pulse  from  its  initial  value 
to  zero;  the  rarefaction  wave  proceeding  forward  is  already  dropping  the 
pressure  such  that  when  the  two  rarefactions  cross  the  pressure  falls  below 
zero.  The  negative  pressure  or  tension  will  be  related  to  the  amount  the 
pressure  has  dropped  behind  the  initial  pulse  when  the  rarefaction  from  the 
front  surface  has  reached  a  given  position.  In  figure  1  the  tension  will  be 
approximately  the  pressure  drop  "h"  that  has  entered  the  material  before  the 
receding  rarefaction  has  reached  the  metal-HE  interface  and  a  pressure  pulse 
has  been  transmitted  back. 

If  either  the  magnitude  or  rate  of  increase  of  tension  is  great  enough,  the 
plate  may  fracture  or  spall. 

THE  CALCULATIONAL  MODEL. 

The  one-dimensional  hydrodynamic  equations  were  solved  by  the  finite 

1 

difference  technique  of  Von  Neumann  and  Richtmeyer.  The  tension  profiles 
for  three  different  plate  thicknesses  were  calculated  using  a  high-speed  com¬ 
puter  together  with  an  assumed  equation  of  state  and  a  method  of  simulating 
the  burning  of  HE.  The  equation  of  state  of  the  metal  in  compression  was 
obtained  from  experimentally  measured  Hugoniot  points.  However,  in  the 
tension  region,  a  Hooke's  law  of  the  form  P  ■  ”1)  was  assumed.  Here 

p  is  the  reference  density  expressed  in  gms  per  cc.  A  perfect  gas  law 

®  2 
equation  of  state  was  used  for  the  HE. 

Since  the  finite  difference  approximation  is  based  on  a  Lagrangian  de¬ 
scription,  the  physical  quantities  of  each  mass  element  are  known  at  any 
instant  of  time.  Consequently  it  is  possible  to  allow  the  material  to  split  at 
a  zone  boundary  by  introducing  the  appropriate  free  surface  boundary 
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conditions.  By  this  technique  material  separation  or  spalling  may  be 
achieved  when  the  tension  falls  below  a  specified  amount.  Other  spall  cri¬ 
teria  based  on  strain  rate,  internal  energy,  momentum  transferred,  etc.  , 
may  be  easily  introduced. 


EXPERIMENTAL  TECHNIQUE. 


To  measure  the  front  surface  velocity  of  HE-driven  plates,  the  usual  pin 
technique  and  optical  methods  were  used.  The  experiments  used  10cm  of 
Composition  B  ignited  by  a  plane  wave  generator.  The  other  dimensions 
were  chosen  to  maintain  a  one-dimensional  model  during  the  time  of  interest. 


To  measure  the  change  in  velocity 
magnetic  method  was  developed.  The 
through  parallel  magnetic  field  lines, 
of  the  plate  is  described  by 

f  -  “7x15 

Xj 


of  the  front  surface  of  a  metal  plate,  a 
plate  is  the  conductor  advancing 
The  voltage  generated  per  unit  length 

V  ■  voltage 

L  s  length  of  plate 

V  s  velocity 

B  *  magnetic  field 


In  order  to  measure  the  voltage  and  hence  the  plate  velocity,  two  small 
wires  are  placed  normal  to  the  plate.  As  the  pltfe  advances  into  the  wires, 
contact  is  maintained  and  the  voltage  as  a  function  of  time  is  recorded  on 
oscilloscopes. 


DISCUSSION. 

By  examining  a  plot  of  the  front  surface  velocity  time  history  of  the  plate, 
successive  accelerations  may  be  discerned.  The  time  interval  between  con¬ 
secutive  pushes  is  equal  to  the  transit  time  of  a  receding  rarefaction  wave 
plus  the  time  of  a  compression  wave.  From  the  previously  determined 
equation  of  state  of  the  metal,  the  shock  velocity  at  the  pressure  given  by  the 
HE  can  be  calculated.  The  velocity  of  the  rarefaction  wave  is  very  nearly 
equal  to  the  shock  velocity,  so  an  easy  calculation  gives  the  time  interval  be¬ 
tween  consecutive  accelerations  of  the  plate  front  surface. 

The  time  interval  calculated  from  this  simple  model  agrees  with  the 
solution  given  by  the  difference  equations  and  also  with  the  experimental 
results  for  thin  metal  plates.  As  the  thickness  of  the  plate  increases,  the 
time  between  consecutive  accelerations  is  longer  than  the  calculated  times. 
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It  is  assumed  that  the  plate  has  spalled  or  fractured  and  the  previous  simple 
calculation  is  no  longer  valid.  A  complete  analysis  of  the  wave  interactions 
together  with  a  criterion  for  fracture  verifies  this  assumption.  The  results 
of  the  finite  difference  calculations  which  take  into  account  the  spall  phenome¬ 
non  correctly  predict  the  time  intervals  between  accelerations  found  experi¬ 
mentally. 

Intuitively,  increasing  the  thickness  of  the  plate  should  have  the  following 
effects: 

1.  A  larger  part  of  the  rarefaction  wave  could  enter  the  material. 

2.  The  pulse  shape  is  changed,  causing  higher  strain  rate's. 

3.  The  time  which  a  specified  plane  in  the  plate  is  under  tension  will 
increase  because  of  the  increase  in  transit  time. 

Any  of  the  above  mechanisms  may  cause  the  maximum  dynamic  strain  in 
the  plate  to  be  exceeded,  thus  allowing  the  material  to  fracture  or  spall. 

Other  forms  of  the  equation  of  state  in  the  negative  pressure  region  as 
well  as  different  spall  criteria  have  been  investigated  on  the  computer.  The 
simple  model  presented  here  illustrated  the  effectiveness  of  the  computer  as 
a  method  of  analysis. 

RESULTS. 

Figures  1.  lA;  2,  2A,  and  3.  3A  illustrate  the  results  of  finite  difference 
calculations  which  do  not  contain  a  spall  criterion  and  the  corresponding 
experimental  results.  Figures  1  and  lA  show  the  results  of  10  cm  of  HE 
driving  a  0.  375  cm  thick  uranium  plate.  The  agreement  between  the  experi¬ 
ment  and  calculation  is  good  and  it  is  concluded  that  the  plate  did  not  spall. 
Figures  2,  2A,  and  3,  3A  show  the  corresponding  results  when  the  plate  has 
a  thickness  of  0.  545  cm  and  0.75  cm.  For  these  cases  the  acceleration  times 
as  measured  by  experiment  do  not  agree  with  the  calculated  times.  The  con¬ 
clusion  is  that  the  material  has  spalled  and  momentum  can  no  longer  be  trans¬ 
ferred.  Thus  the  spalled  portion  of  the  plate  flies  off  with  the  momentum 
characteristic  of  the  pressure  from  the  head  of  the  original  shock  minus  the 
amount  of  momentum  transferred  back  by  tension  before  spall  occurred. 

The  decrease  in  front  surface  velocity  between  accelerations  shown  by 
the  calculations  results  from  the  tensions  in  the  plate.  Experimentally  this 
is  not  seen,  since  only  average  velocities  are  measured  in  the  time  interval 
between  accelerations. 
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The  calculations  were  repeated  using  a  spall  criterion  of  -0.  50  x  10^^ 

dynci  ;  i.  e.  ,  the  material  was  allowed  to  separate  if  the  tension  at  any  point 
cmr 

in  the  material  exceeded  this  amount.  Figures  4  and  5  show  the  results  of 
these  calculations  as  compared  to  experiment. 

Similar  calculations  using  other  thicknesses  of  plates,  and  a  spall  criterion 
based  upon  tension  of  about  050  x  10^^  cm^^  agreed  with  the  experiment. 

The  assumed  equation  of  state  in  the  region  of  negative  pressure  and  the  cri¬ 
teria  for  spall  used  are  certainly  not  correct  in  every  detail;  however,  the 
overall  good  agreement  with  experiment  is  justification  for  applying  the  tech¬ 
nique  to  other  HE  metal  plate  systems. 

CONCLUSION. 

As  stated  in  the  first  section,  spall  results  when  two  rarefactions  cross. 

An  obvious  method  to  defeat  spall  would  be  to  cancel  one  of  the  rarefactions. 

To  accomplish  this  a  void  can  be  placed  between  the  HE  and  the  metal.  The 
transmitted  shock  into  the  material  now  has  a  rising  profile  instead  of  a  pro¬ 
file  that  falls  behind  the  shock  front.  Figures  6  and  6A  show  the  results  of 
calculation  and  experiment  where  a  1.  5  cm  void  was  introduced  between  the 
0.  75  cm  plate  and  the  HE. 

The  experimental  methods,  while  accurate  for  measuring  the  changes  in 
front  surface  velocity,  only  give  average  velocities  during  the  intervals  be¬ 
tween  accelerations.  It  would  take  separate  experiments  to  establish  the 
slowing  down  of  the  front  surface  as  predicted  by  the  calculations.  Experi¬ 
ments  of  this  type  are  presently  being  conducted. 


(This  document  was  originally  released 
as  UCRL  6356) 
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Figure  4  Figure  4-A 
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DYNAMIC  LARGE  DEFLECTION  OF  SHELL  STRUCTURES 

Dr.  W.  Herrmann 
T.  H.  H.  Plan 
J.  W.  Leech 

Aeroelastic  and  Structures  Research  Laboratory 
Massachusetts  Institute  of  Technology 

The  problem  considered  concerns  large  deformations  of  shell-type 
structures  subjected  to  non-uniformly  distributed  intense  loads  of  very  short 
duration.  A  brief  review  will  be  given  of  the  methods  and  assumptions  on 
existing  dynamic  analysis  of  elastic  and  plastic  structures. 

Where  the  impulse  is  sufficiently  low,  only  elastic  deformations  result. 
The  problem  of  lateral  motion  (normal  to  the  shell  surface)  is  somewhat  more 
difficult  than  that  of  the  longitudinal  motion.  In  the  case  of  elastic  deforma¬ 
tions.  the  equation  of  lateral  motion  involves  a  fourth-order  differential  equa¬ 
tion  while  the  equation  for  longitudinal  motion  involves  only  a  second-order 
differential  equation. 

Transient  solutions  of  lateral  elastic  motion  can  be  obtained  hy  the 
classical  normal  mode  solution  or  by  the  wave  solution.  When  the  duration 
of  the  applied  dynamic  load  is  not  too  short,  the  normal  mode  method  is  the 
suitable  choice.  However,  for  problems  which  involve  discrete  pulse  loads 
of  very  short  duration,  the  convergence  of  the  normal  mode  solution  becomes 
questionable  and  the  wave  solution  must  be  adopted.  It  has  been  shown,  how¬ 
ever,  that  if  elementary  beam  theory  is  used  the  resulting  flexual  wave 
propagation  leads  to  dispersion.  The  resulting  phase  velocity  is  found  to  be 
inversely  proportional  to  the  wavelength.  For  a  wave  of  infinitely  short 
wavelength  the  propagation  velocity  is  infinite.  One  can  thus  conclude  that 
elementary  beam  theory  is  inadequate  for  short  impulsive  loads  because  it 
leads  to  the  physically  impossible  conclusion  that  disturbances  are  propagated 
instantaneously,  throughout  the  beam.  Refinements  to  elementary  theory 
have  been  made  by  Rayleigh  and  by  Timoshenko  by  taking  into  account  the 
rotating  inertia  of  tbe  beam  cross  section  and  the  transverse  shear  deforma¬ 
tion.  From  this  improved  theory  it  can<be  shown  that  there  are  two  types  of 
waves  propagated  along  the  beam,  one  of  which  involves  bending  moment  and 
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cross-section  rotation  while  the  other  involves  transverse  shear  force  and 
lateral  displacement.  The  two  groups  are  traveling  at  different  propagation 
velocities  but  they  are  coupled  with  each  other,  and  the  shape  of  each  wave 
changes  completely.  The  coupled  wave  equations  can  be  solved  by  the  method 
of  characteristics. 

Studies  of  plastic  deformation  of  beams  caused  by  lateral  impulsive  load¬ 
ing  are  of  recent  origin,  having  been  initiated  in  this  decade.  Most  of  the 
studies  have  been  based  on  simple  beam  theory  with  the  moment-curvature 
relationship  being  considered  as  elastic-plastic  (figure  l.a)  or  rigid- plastic, 
(figure  l.b)  and  the  plasticity  being  described  generally  as  either  perfectly 
plastic  (figure  1.  a  and  b)  or  as  linear  strain-hardening  (figure  1.  c). 

Analyses  of  elastic-plastic  beams  have  also  been  made  by  using  the 
normal  mode  approach.  In  this  analysis  the  solution  must  be  divided  into 
individual  phases.  When  the  bending  moment  at  any  point  reaches  the  yield 
bending  moment,  M^,  a  single  plastic  hinge  is  formed  and  in  the  subsequent 
period  the  beam  must  be  considered  as  two  elastic  beams  connected  by  the 
hinge.  By  considering  the  proper  continuity- relations  between  displacements 
and  velocities  at  the  joint,  the  solution  can  be  established. 

Solutions  of  lateral  motion  of  plastic  beams  have  also  been  obtained  using 
the  wave  propagation  approach.  The  solution  by  Duwez,  Clark,  and  Bohnenblust 
is  based  on  simple  beam  theory  while  the  solution  by  Plass  is  based  on  an  im¬ 
proved  theory  taking  into  account  the  rotational  inertia  and  the  shear  deforma¬ 
tion.  These  authors,  however,  consider  only  the  problem  of  a  semi-infinite 
beam  subjected  to  a  lateral  disturbance  at  one  end. 

All  the  work  on  plastic  beams  to  date  has  emphasized  the  permanent 
deformation  due  to  the  impulsive  loading.  For  the  cases  where  the  plastic 
deformation  is  much  larger  than  the  elastic  deformation,  the  simplified 
rigid-plastic  analysis  is  justified.  From  the  basic  assumption  of  negligible 
elastic  strain  it  can  be  concluded  that  the  motion  of  the  beam  involves  only 
rigid  body  motions  of  various  beam  segments  which  are  connected  by  plastic 
hinges.  This  type  of  analysis  was  first  developed  by  Lee  and  Symonds  and 
was  followed  by  other  authors.  A  simple  example  using  rigid- plastic  analysis 
is  given  in  the  following. 

Consider  a  simply  supported  beam  acted  on  by  a  uniformly  distributed 
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impulsive  load  p{t).  (Figure  2)  When  the  applied  load  is  small,  such  that 
the  maximum  bending  m  jment  is  less  than  the  yield  bending  moment,  the  beam 
remains  rigid.  The  maximum  bending  moment  which  occurs  at  the  mid-span 
is  equal  to  .  Thus  when  p  •  8  M  ,  a  plastic  hinge  is  formed  at  the  mid- 

span.  If  the  applied  load  further  increases,  the  beam  will  move  with  two  rigid 
segments  rotating  about  the  end  supports.  The  magnitude  of  the  angular 
acceleration  of  the  rigid  segments  must  be  such  that  the  bending  moment  at  the 
mid- span  is  maintained  at  .  From  figure  2  we  obtain 

^  -  M 

8  24  ° 

or  'ft  =  -  M  ) 

m  X-  8  o 

When  the  impulsive  load  is  increased  continuously,  a  point  will  be  reached 
when  the  location  of  maximum  bending  moment  begins  to  spread  outward  from 
the  mid- span  point.  This  is  the  third  phase  of  the  motion,  when  two  plastic 
hinges  begin  to  appear  and  the  beam  will  move  with  three  rigid  segments.  Two 
of  these  segments  rotate  about  the  end  support  while  the  segment  at  the  center 
undergoes  uniform  translation  (figure  2.  f.  )  The  plastic  hinges  are  moving 
away  from  the  center  when  the  load  is  increasing  and  they  are  moving  toward 
the  center  when  the  load  has  been  released.  The  rigid- plastic  analysis  involves 
essentially  the  determination  of  the  motion  in  the  various  phases  until  the  ve¬ 
locity  of  the  beam  becomes  zero  everywhere,  and  the  final  plastic  deformation 
can  be  evaluated.  The  propagation  of  the  plastic  hinge  in  the  rigid  plastic 
analysis  is  not  of  the  same  nature  as  the  propagation  of  elastic  or  plastic  waves. 
In  the  present  case,  there  is.  in  fact,  no  inertia  associated  with  the  moving 
plastic  hinge. 

Another  refinement  of  the  plastic  beam  analysis  is  the  inclusion  of  the 
axial  load  effect  on  the  yield  bending  moment.  For  a  beam  restrained  at  both 
ends  from  longitudinal  displacements,  axial  tensile  stress  is  developed  when 
the  deflection  of  the  beam  becomes  larger  than  the  thickness.  For  a  rigid- 
plastic  beam  of  rectangular  cross  section  the  magnitudes  of  the  bending  moment , 
M.  and  the  axial  force,  N  .  at  the  yield  condition  are  related  by  the  following 
interaction  equation: 
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where  s  yielding  bending  moment  when  the  axial  load  N  is  zero  and 
a  yielding  axial  load  when  the  bending  moment  M  is  zero. 

It  can  be  seen  that  in  the  initial  phase  where  the  deflection  and  the  axial 
load  are  small,  the  plastic  work  involves  essentially  bending  at  the  plastic 
hinge.  However,  when  the  deflection  increases,  the  axial  tensile  load  also 
increases,  and  the  plastic  work  will  involve  mainly  axial  strain,  instead,  of 
bending  strain.  In  the  limiting  case  where  the  axial  load  becomes  equal  to 
,  the  beam  becomes  a  plastic  string  under  constant  tension  .  (Thus 
this  last  phase  must  be  treated  differently).  Such  an  analysis  was  made  by 
Symonds  and  Mentel.  A  treatment  of  a  low  arch  under  lateral  load  was  made 
by  Chen,  Hsu  and  Plan.  In  the  latter  case  axial  compressive  stress  is  de¬ 
veloped  when  the  arch  deforms  inward. 

All  of  the  previously  described  plastic  theories  apply  to  small  deflections 
in  that  geometrical  nonlinearities  are  excluded;  and  in  addition,  they  make 
further  simplifying  assumptions,  either  neglecting  elastic  effects  by  assuming 
rigid-plastic  behaviour,  or  considering  only  infinitesimal  deflections  for  which 
the  axial  (membrane)  stress  is  negligible. 

Extension  to  large  deformations  involves  not  only  the  introduction  of 
axial  (membrane)  stresses  but  additional  nonlinearities  caused  by  the  geo¬ 
metrical  effects.  The  resultant  equations  require  numerical  solution.  It  was 
decided  to  retain  elastic  effects  in  this  analysis  since  certain  evidence  was 
found  to  suggest  that  the  energy  stored  in  elastic  strain  in  certain  deformed 
shell  configurations  is  not  negligible.  In  particular  an  appreciable  amount  of 
energy  appears  to  be  stored  in  elastic  strain  in  circular  cylindrical  and 
spherical  shells  when  these  have  been  deformed  to  the  order  of  half  the  radius, 
as  may  be  shown  by  cutting  the  shell  at  some  convenient  location  and  measur¬ 
ing  the  work  required  to  close  the  gap. 

The  method  of  solution  requires  setting  the  equation  of  motion  for  the 
shell,  written  for  large  deflections,  together  with  the  stress-strain  relations, 
in  finite  difference  form.  The  resulting  initial  value  problem  is  solved, 
subject  to  the  boundary  conditions  describing  the  constraints  on  the  shell.  The 
load  may  be  introduced  as  a  varying  boundary  condition  if  the  loading  is  of 
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relatively  long  duration.  Alternatively  if  the  load  duration  is  very  short,  and 
the  structure  does  not  deform  appreciably  during  the  time  of  load  application, 
the  load  may  be  characterised  as  an  impulse,  giving  rise  to  an  initial  velocity 
distribution  in  the  shell. 

Initially  only  cylindrical  shells  (rings)  are  being  considered,  thus  elimin¬ 
ating  one  space  dimension.  The  ring  need  not  be  initially  circular.  The  shell 
may  also  be  open  ended,  with  appropriate  end  constraints,  and  therefore 
equivalent  to  a  curved  beam. 

The  equilibrium  equation  in  finite  difference  form  may  be  described  in 
physical  terms  as  follows.  The  ring  is  considered  to  be  composed  of  a 
number  of  links  joined  by  hinges.  The  links  are  assumed  rigid  in  bending  but 
elastic-plastic  in  extension.  The  hinges  are  considered  to  be  inextensible, 
but  elastic-plastic  in  bending.  The  mass  of  the  ring  is  divided  into  discrete 
masses  concentrated  at  the  hinges.  (Figure  3.  ) 

As  in  the  finite-deflection  analysis  of  Symonds  and  Mentel,  mentioned 
previously,  there  is  an  Interaction  between  bending  moment  and  axial  force 
at  the  yield  condition.  For  the  elastic- perfectly  plastic  stress- strain  relation 
used  here,  the  interaction  equation  cannot  be  written  simply,  for  a  rec¬ 
tangular  cross  section.  Particular  difficulty  is  experienced  when  stress 
reversals  occur.  The  difficulties  are  circumvented  by  idealising  the  cross 
section  into  an  I  beam  (or  sandwich  plate)  of  equivalent  force  and  moment 
carrying  capacity.  The  flanges  are  considered  to  behave  in  an  elastic- 
perfectly  plastic  manner,  and  a  relatively  simple  relation  between  extensional 
strain  and  curvature,  and  axial  stress  and  moment,  can  be  found. 

Since  stress  reversals  may  occur,  and  unloading  will  occur  elastically,  a 
multivalued  relation  results,  the  correct  relation  depending  on  the  previous 
loading  history.  It  is  therefore  necessary  to  provide  the  correct  logic  in  the 
program  to  keep  track  of  the  loading  history,  in  order  to  choose  the  correct 
relation  between  the  axial  stress  and  strain,  and  the  moment  and  curvature. 

Preliminary  results  for  the  case  of  a  beam  subjected  to  a  uniform  initial 
velocity  indicate  that  traveling  hinges  are  formed,  as  in  the  rigid- plastic 
case  of  Lee  and  Symonds.  It  is  anticipated  that  the  program  will  give  insight 
into  the  failure  modes  to  be  expected  in  various  configurations,  and  information 
about  the  relative  amounts  of  energy  absorbed  in  elastic  and  plastic  deformation, 
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MOTION  OF  AN  ELASTIC  HALF-SPACE 

Dr.  C,  M.  Ablow 

Dr.  Roy  C.  Alverson 
Stanford  Research  Institute 

Underground  shelters  are  designed  to  withstand  the  motion  of  the  earth 
under  blast  loadings  caused  by  explosions  on  or  above  the  earth's  surface. 

A  first  approximation  to  that  motion  is  obtained  by  computing  the  small 
motions  of  an  isotropic  elastic  half  space  under  circularly  symmetric  applied 
tractions.  Features  of  the  motion  of  major  importance  are  preserved  by  a 
computing  method  which  accurately  predicts  the  strength  and  time  of  arrival 
of  the  first  seismic  shock.  Such  a  method  is  at  hand  when  shear  and  com- 
pressional  wave  motions  are  separated  and  independent  variables  are  constant 
on  wave  fronts,  i.  e.  when  characteristic  coordinates  are  introduced. 

For  axisymmetric  motions  of  the  elastic  half  space  a  cylindrical  co¬ 
ordinate  system  (r,  T.z)  is  appropriate.  It  is  assumed  that  there  is  no  ro¬ 
tation  about  the  axis  so  that  all  functions  appearing  are  independent  of  f.  It 
is  then  known^  that  the  equations  of  motion  and  of  Hook's  lew  may  be  reduced 
to  the  solution  of  two  equations  for  two  potential  functions  cp  and  P,  function 
0  determining  the  compression  and  function  0  the  shear: 

,  2 

(P  +(p/r  +  ffl  s  V1.1./0 
rr  ^r'  ^zz  ^tt' 

e  +e/r+e  -e/r^ 

rr  r'  zz  '  tt' 

where  a^«{X.+2n)/p,  \  and  |i  are  the  Lame  elastic  parameters, 

and  p  is  the  density.  Although  the  reduction  to  the  above  form  is  carried 
through  in  the  reference  only  for  constant  (x,  p,  the  same  manipulations 
are  permissible  if  \  varies  as  a  function  of  z.  Thus  the  present  analysis 
holds  for  an  earth  of  constant  density  and  shear  wave  speed  but  compression- 
al  wave  speed  varying  with  depth.  Displacements  in  the  radial  direction,  u, 
and  in  the  axial  direction,  w,  are  determined  by  ro  and  0  to  be 

u«(p  -0  ,  w»cp  +  0  +  0/  r. 

r  z  XT' 

1  W.  M.  Ewing,  W.  S.  Jardetzky,  and  F.  Press,  Elastic  Waves  in  Layered 
Media,  McGraw  Hill,  New  York,  1957,  Chapter  I 
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If  the  surfaces  in  (r,  z,  t)  space  on  which  function  f  is  constant  are  shear 
wave  fronts  then  f  is  a  solution  of  the  characteristic  equation 

-  P'  (f/  ^  i^)  . 

Fundamental  shear  wave  front  surfaces  are  the  right  circular  cones  with 
arbitrary  vertex,  axis  parallel  to  the  time  axis,  and  half  vertex  angle  arc 
ctn  p .  It  may  be  shown  that  any  shear  wave  front  is  an  envelope  of  such 
cones.  For  the  characteristic  coordinate  surface  c  >  0  there  was  taken  the 
envelope  of  cones  with  vertices  on  the  line  marking  the  leading  edge  of  the 
disturbance  on  z  ■  0,  the  surface  of  the  earth.  The  location  of  c  -  0  in 
(r,  z,  t)  space  marks  the  progress  of  the  initial  shear  disturbance  throughout 
the  medium. 

The  two  other  characteristic  coordinate  planes  a  b  0  and  S  =  0  were 
taken  as  the  surfaces  representing  plane  waves  parallel  to  the  surface  of  the 
earth,  being  at  that  surface  at  time  zero,  and,  respectively,  going  out  into 
the  earth  or  returning  from  the  depths.  Thus  a  a  0  is  one  planar  envelope  of 
the  cones  with  vertices  on  the  r-axis  in  (r,  z,  t)  space  while  S  s  0  is  their 
other  planar  envelope. 

The  same  definitions  give  the  a  s  0,  b  ■  0,  and  c  a  0  characteristic  sur¬ 
faces  for  the  compressional  wave  equation,  except  that  the  variability  of  a 
with  depth  distorts  the  right  circular  cones  into  surfaces  of  somewhat  similar 
shape  called  conoids. 

To  obtain  a  complete  coordinate  system  in  (r,  z,  t)  space,  the  surface 
c  ■  k  was  taken  as  the  c  ■  0  surface  moved  k  units  parallel  to  itself  in  the 
direction  of  the  t-axis,  and  similarly  for  the  other  coordinates.  Analytically, 
if  f(r,  z,  t)  is  that  solution  of  the  characteristic  equation  for  which  f(r,  z,  t)  a 
0  represents  the  c  >  0  surface,  then  the  equation  for  c  b  k  was  taken  to  be 
f(r,  z,  t-k)  s  0. 

It  is  evident  that  there  is  a  high  degree  of  arbitrariness  in  the  choice  of 
characteristic  coordinate  system.  The  choice  made  above  makes  coordinate 
planes  of  the  initial  compressional  and  shear  wave  fronts  and  also  provides 
simple  plane  surfaces  in  characteristic  space  for  the  surfaces  r  s  0  and 
z  s  0  on  which  boundary  conditions  are  applied  as  well  as  the  surfaces 
t  a  constant  on  which  the  computed  output  is  conveniently  displayed.  Further, 
on  the  boundary  z  >  0,  where  compressional  and  shear  waves  satisfy  combined 
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conditions,  one  has  a  s  a,  b  :  S,  and  c  ■  c  so  that  there  is  no  interpolation 
problem. 


Physical  surface 

Characteristic  surface  ’ 

for  compression 

for  shear 

Initial  wave 

o 

II 

u 

C  s  0 

r  =  0 

a  s  c 

a  X  c 

Z  s  0 

a  s  b 

II 

t  .  t^ 

a  +  b  s  2tj 

a  +6  ■  2tj 

Setting 


<Pa  ■ 

■  B, 

s  C. 

e;.  A. 

9^  s  B,  and 

6-  .  C 
c 

permits  the  second>order  potential  equations  to  be  replaced  by  the  first' 
order  systems 

4A.+2(ac  +  DA  -2(ac  -  1)B  -  a^V^cC  .  0 

D  Z  C  Z  C 

A.  ■  B  ,  B  ■  C.  I 

b  a  c  b 

and 

4A  c  +  2(pc  +  DA-  -  2(pc  -  DBp  -  +  p^e/r^  .  0  , 

D  Z  C  Z  C 

A£  -  Bj  .  B-  ,C^ 

where  the  first  set  of  equations  replaces  the  potential  equation  for  the  di- 

latational  waves  and  the  second  set  replaces  the  potential  equation  for  the 

shear  waves.  In  these  equations  V^c  stands  for  c  +  c  /r  +  c  ,  and  it, 

rr  r  zz 

as  well  as  the  other  quantities  appearing,  needs  to  be  expressed  in  the  char¬ 
acteristic  variables  before  the  integrations  can  begin. 

Finite  difference  approximations  to  the  differential  equations  have  been 
written  using  a  cubical  lattice  in  each  characteristic  coordinate  system.  As 
is  to  be  expected  with  a  hyperbolic  system  in  characteristic  form,  the  equa¬ 
tions  are  explicit,  the  solution  is  obtained  by  a  process  of  marching  through 
the  nodes  of  the  lattice,  and  stability  of  this  process  is  assured. 
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The  initial  wave  fronts,  being  coordinate  surfaces,  are  sharply  defined 
in  the  finite  difference  calculation.  No  jump  condition  is  needed  across  the 
wave  front  as  the  proper  discontinuity  is  preserved  in  the  finite  difference 
forms. 

The  finite  difference  computation  will  smear  sharp  loading  or  unloading 
waves  behind  the  initial  front  unless  they  fall  along  a  characteristic  co¬ 
ordinate  surface,  i.  e.  ,  are  plane  waves  parallel  to  the  surface  of  the  earth 
or  are  waves  of  the  same  speed  and  extent  as  the  initial  wave  delayed  in  time. 

The  main  disadvantage  of  the  characteristic  method  is  that  three  co¬ 
ordinate  systems  need  to  be  carried  along  together,  (r,  z,  t),  (a,  b,  c),  and 
(a,  6,  c). 

Displacements,  velocities,  and  stresses  in  three  different  models  of  the 
earth  under  airburst  loading  are  being  computed  and  will  appear  in  forth¬ 
coming  reports.  The  three  models  all  consider  the  earth  as  an  isotropic, 
elastic,  half  space.  In  the  first  model  the  half  space  is  homogeneous,  in 
the  second  it  is  homogeneous  with  a  homogeneous  overburden,  and  in  the 
third  the  compressional  wave  speed  vhries  linearly  with  depth. 
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GENERAL  COMMENTS  ON 
NUMERICAL  INTEGRATION  SCHEMES 

Dr.  John  G.  Trullo 
Lawrence  Radiation  Laboratory 

With  regard  to  the  numerical  integration  errors  inherent  in  the  von 
Neumann-Richtmyer  equations,  and  the  effect  of  these  errors  in  spall  calcula¬ 
tions  *  I  would  like  to  make  the  following  comments. 

First,  as  stated,  the  linear  Q  has  the  effect  of  eliminating  oscillations  be¬ 
hind  a  shock.  It  introduces  much  more  damping  than  the  quadratic  Q.  When 
velocity  differences  from  zone  to  zone  are  not  large,  squaring  them  tends  to 
produce  Q's  so  small  as  to  be  ineffective  for  damping  oscillations  in  density, 
etc. 

However,  a  basic  problem  in  describing  rarefactions  is  that  any  signal 
which  does  not  move  with  the  velocity  of  the  grid  will  diffuse.  This  appears  in 
the  rounding  of  the  head  of  a  rarefaction  wave,  and  throws  the  timing  off  in 
rarefaction  interaction  problems.  For  most  problems  this  is  not  a  very  im¬ 
portant  effect.  However,  it  is  important  in  any  code  calculation  of  the  depth 
of  spall,  which  is  one  of  the  main  Jobs  of  the  (Boeing)  code. 

If  a  shock  arrives  at  a  free  surface,  the  head  of  the  reflected  rarefaction 
wave  may  be  6  to  8  zones  deeper  in  the  material  than  it  should  be.  (The 
error  in  the  position  of  the  rarefaction  head  is,  of  course,  characterized  by 
a  fixed  number  of  zones,  and  not  a  fixed  distance. )  An  error  of  6  or  8  zones 
is  not  serious  if  there  are  600  zones  in  the  spalled  piece,  but  there  is  ordin¬ 
arily  a  limit  to  the  number  of  zones  in  these  problems  —  a  practical  limit  in 
calculation  time.  Also,  a  method  which  requires  on  the  order  of  1,  000  zones 
to  produce  an  acceptable  solution  to  a  one -dimensional  problem  is  perhaps  not 
the  right  method  to  use  on  the  problem.  For  that  reason,  other  methods  are 
being  investigated. 

The  various  numerical  integration  procedures  have  application  to  a  series 
of  experiments  which  have  been,  and  are  being,  carried  out  to  investigate  the 

"Butler,  Gunning,  Jr.,  and  Young,  Daniel  M. ,  On  the  Treatment  of  Rarefaction 
Using  a  Dissipative  Hydrodynamics  Code 

Wilkins,  Mark,  Calculation  of  Spall  Based  on  a  One -Dimensional  Model 
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mechaniim  of  the  spall  process.  The  purpose  of  the  experiments  is  to  de¬ 
termine  whether  or  not  the  spall  mechanism  is  time  dependent.  The  experi¬ 
ments  are  based  on  the  scaling  laws  of  hydrodynamics,  as  described  in  the 
proceedings  of  last  year's  meetings  (AFSWC-TR-60-12).  In  each  experiment, 
one  plate  is  driven  at  another  plate.  The  absolute  plate  thicknepses  are  varied 
from  experiment  to  experiment,  keeping  constant  the  ratio  of  plate  thicknesses 
and  the  relative  velocity  of  the  plates  before  impact.  Then  the  pressure  dis¬ 
tribution  in  the  plates  is  the  same  in  each  experiment  at  corresponding  positions 
md  times.  Therefore,  if  there  is  no  delay  time  to  fracture,  one  should  observe 
threshold  spall  for  the  same  relative  velocity  of  the  plates,  regardless  of  their 
absolute  thicknesses.  So  far,  this  does  not  seem  to  be  the  case,  although  we 
have  drawn  no  firm  conclusion  as  yet.  Several  sets  of  experiments  have  been 
done,  for  three  or  four  plate  thickness  ratios.  Results  have  been  obtained  both 
by  exploding  foil  and  air  gun  techniques;  they  agree  rather  closely. 

If  the  spall  mechanism  is  time  dependent,  then,  in  order  to  calculate  the 
depth  of  spall,  etc. ,  as  part  of  the  problem  of  motion,  it  will  be  necessary  to 
build  into  the  code  the  experimental  delay  time  to  fracture,  as  a  function  of  the 
applied  stress.  For  stresses  and  strains  outside  the  linear  elastic  range,  the 
deduction  of  this  delay  time  from  the  experimental  data  leads  to  the  further 
problem  of  determining  the  stress-strain  relation  for  large  tensile  stresses. 

If  spall  occurs  instantly  at  a  fixed  teniile  stress,  the  appropriate  code  modifi¬ 
cation  is  much  simpler.  In  either  case,  a  quantitative  empirical  (or  semi- 
empirical)  account  of  one -dimensional  spall  will  result  from  experiements  of 
the  kind  described.  A  technique  for  integrating  the  equations  of  motion  is 
therefore  acceptable  only  if  it  allows  us  to  make  full  use  of  this  knowledge  to 
predict  the  occurrence  of  spall,  the  depth  of  spall,  and  other  details  of  the 
spall  process.  On  this  basis,  the  equations  of  von  Neumann  and  Richtmyer, 
unmodified,  are  inadequate. 

I  would  like  to  comment  also  on  an  accuracy  limitation  inherent  in  any 
numerical  scheme  which  integrates  the  hydrodynamic  equations  by  stepwise 
advances  in  time.  Many  attempts  have  been  made  to  increase  the  accuracy  of 
the  equations  by  "high-order  differencing.  "  That  is,  the  difference  analogs 
of  the  hydrodynamic  equations  are  written  in  such  a  way  that  the  error  in¬ 
troduced  in  a  dependent  variable  in  a  single  timestep  is  of  high  order  in  the 
sense  of  a  Taylor's  Series  expansion  in  powers  of  the  space-  and  time-steps. 


182 


TN-61-29 


Efforts  in  this  direction  have  had  very  limited  success  to  my  knowledge.  There 
is,  I  think,  a  simple  and  basic  reason  for  this:  quite  apart  from  shocks  and 
contact  discontinuities,  the  solutions  of  the  equations  of  motion  are  not  analytic 
for  flows  of  even  moderate  complexity.  For  example,  in  one-dimensional 
motions  of  practical  interest,  the  x  -  t  plane  is  criss-crossed  by  curves  across 
which  the  first  derivatives  of  the  sound  speed  or  particle  velocity  with  respect 
to  X  or  t  are  discontinuous.  Then  the  introduction  of  high-order  differences 
into  the  difference  equations  can  actually  bring  about  a  decrease  in  this  accuracy 
of  the  equations,  because  the  derivatives  to  which  these  differences  tend  in  the 
limit  of  fine  zoning  are  infinite  at  various  points  in  the  flow.  Although  the 
truncation  error  in  a  single  time-step  is  reduced  at  most  points  by  high-order 
differencing,  the  error  introduced  at  the  relatively  few  points  on  the  singular 
lines  of  the  flow  can  more  than  offset  the  gain  in  accuracy  elsewhere.  It  seems 
that  second-order  accuracy  is  about  the  most  that  can  or  should  be  asked  of  a 
general  purpose  hydrodynamic  code  which  operates  through  stepwise  advances 
in  time.  To  get  around  this  difficulty,  and  thereby  write  rapidly  convergent 
difference  equations,  one  must  see  to  it  that  differences  are  never  taken  across 
singular  surfaces  of  a  flow  in  the  space -time  continuum.  Since  the  singular 
surfaces  are  characteristic  surfaces  (apart  from  shocks  and  contact  discon¬ 
tinuities),  the  difficulty  is  overcome  by  using  a  set  of  characteristic  coor¬ 
dinates  as  independent  variables.  This,  it  seems  to  me,  is  the  basic  reason  for 
the  much  more  rapid  convergence  obtained  by  differencing  the  equations  in 
characteristic  form.  However,  the  use  of  the  characteristic  equations  for 
numerical  purposes  has  its  own  difficulties,  especially  in  the  treatment  of 
shocks  and  contact  discontinuities,  and  in  two  or  three  space  dimensions. 

Since  the  flow  of  a  compressible  fluid  is  a  problem  in  wave  motion 
(although  nonlinear),  it  is  natural  to  examine  the  question  of  accuracy  from  the 
point  of  view  of  a  Fourier  expansion  as  well  as  a  Taylor's  Series.  Along  these 
lines,  the  simple  problem  of  linear  wave  propagation  relative  to  a  coordinate 
mesh  has  been  studied.  For  concreteness,  the  motion  can  be  thought  of  as  the 
uniform  translation  of  a  bar  of  material  with  variable  density  and  zero  pressure, 
through  an  Eulerian  grid.  The  difficulty  is  that  the  density  is  known  basically 
at  the  Eulerian  zone  centers,  while  the  mass  flux  needed  to  calculate  changes 
in  density  must  be  defined  at  the  zone  boundaries.  There  are  many  ways  in 
which  the  mass  flux  terms  can  be  calculated  from  the  basic  densities,  each 
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corresponding  to  a  particular  construction  of  a  complete  density  profile  from 
the  basic  discrete  density  field.  Apart  from  some  trivial  cases,  different 
definitions  of  the  mass  flux  will  produce  different  cycle-by>cycle  density 
changes.  Thus,  instead  of  the  correct  undisturbed  propagation  of  a  density 
profile,  the  density  profile  will  be  distorted  in  the  numerical  calculations.  It 
is  worth  noting  that  this  transport  problem  appears  whenever  a  disturbance 
propagates  at  a  finite  speed  relative  to  the  coordinate  system  used  to  describe 
it.  In  a  Lagrangian  frame,  the  problem  arises  for  sound  wave  propagation. 

One  simple  definition  of  the  density  used  to  calculate  the  mass  flux  at  a 
zone  boundary  is  the  so-called  "backward"  definition,  i.  e. ,  take  the  nearest 
zone-centered  density  in  the  direction  corresponding  to  the  tail  of  the  velocity 
vector  at  the  zone  boundary.  If  this  definition  is  used,  it  is  found  that  the 
amplitude  of  an  8 -zone  sine  wave  decays  by  a  factor  of  10  by  the  time  the  wave 
has  travelled  two  or  three  wavelengths  relative  to  the  grid.  A  cycle-by-cycle 
analysis  of  the  results  shows  clearly  why  this  rapid  diffusion  takes  place;  the 
point  will  not  be  pursued  here.  The  rate  of  diffusion  decreases  as  more  zones 
are  used  per  unit  wavelength,  so  that  the  problem  centers  mainly  on  short 
waves. 

To  do  a  better  job  with  sine-wave  propagation,  another  simple  definition 
of  the  transport  density  was  adopted.  This  is  the  so-called  "centered"  defin¬ 
ition,  i.  e. ,  the  zone-centered  densities  on  either  side  of  the  zone  boundary  are 
averaged  (a  small  correction  term  is  added  for  purposes  of  numerical  stabil¬ 
ity).  With  this  definition,  the  8-zone  wave  propagates  with  very  little  distortion 
and  no  noticeable  diffusion.  However,  if  the  density  profile  is  a  step  function, 
then  the  centered  differencing  scheme  quickly  leads  to  a  train  of  large  oscilla¬ 
tions  behind  the  step.  The  growth  of  these  oscillations  is  easy  to  understand 
by  following  the  calculation  time-step  by  time-step.  The  initial  spurious  dis¬ 
turbance  is  a  zone-by-zone  oscillation.  Thus,  the  problem  again  centers  on 
the  description  of  very  short  wavelengths.  In  a  finite  difference  grid,  the 
shortest  possible  sine  wave  is  spread  over  four  zones.  These  shortest 
veives,  and  even  somewhat  longer  ones,  are  very  coarsely  defined,  and  do  not 
propagate  correctly  through  the  grid.  To  eliminate  errors  from  this  source, 
an  obvious  stratagem  is  to  suppress  selectively  all  sinusoidal  components 
shorter  than  some  predetermined  wavelength,  without  disturbing  the  others. 

One  way  this  has  been  done,  according  to  Dr.  Cecil  Leith  at  Livermore,  is 
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actually  to  Fourier  analyze  the  whole  signal  and  subtract  out  the  unwanted 
components.  We  have  another  scheme  that  accomplishes  essentially  the  same 
thing  -  waves  spread  over  any  particular  numbers  of  zones  can  be  damped  out. 
By  suppressing  the  saw-toothed  waves,  all  oscillations  are  wiped  out  but  the 
first  two,  for  a  uniformly  translating  step  function.  The  amplitude  of  the  first 
oscillation  behind  the  step  is  only  slightly  reduced,  while  the  amplitude  of  the 
second  oscillation  is  about  halved.  By  suppressing  the  waves  whose  half-width 
is  three  zones  or  less,  the  second  oscillation  is  eliminated,  and  the  amplitude 
of  the  first  oscillation  is  reduced  by  about  a  factor  of  5  (relative  to  the 
amplitude  characteristic  of  centered  differencing  with  no  suppression).  At  the 
same  time,  sine  waves  propagate  exactly  as  in  the  case  of  centered  differenc¬ 
ing.  The  method  is  therefore  fairly  satisfactory  for  wave  propagation  in  a 
finite  grid. 

The  wave  propagation  study  was  undertaken  largely  to  provide  a  means  of 
handling  the  transport  terms  which  appear  when  the  equations  of  motion  are 
written  in  a  non-Lagrangian  coordinate  system.  The  treatment  of  transport 
resulting  from  this  study  has  been  put  into  a  code  which  describes  the  one¬ 
dimensional  linear  motion  of  a  single  material  with  two  Lagrangian  boundaries, 
one  of  which  is  not  allowed  to  move.  Between  the  boundaries  at  any  given  time 
are  a  fixed  number  of  zones  of  equal  width.  The  grid  is  therefore  neither 
Lagrangian  nor  Eulerian;  the  zone  boundaries  do  not  move  with  the  particle 
velocity  and  they  are  not  fixed  in  space.  The  most  interesting  problem  run  on 
the  code  to  date  is  that  of  a  steady  shock  and  its  reflection  from  a  rigid  wall. 
The  error  in  the  pressure  behind  the  transmitted  shock  is  about  50  percent 
greater  than  the  error  using  the  von  Neumann-Richtmyer  equations.  This  is 
actually  better  than  it  looks  because  the  mass  of  a  zone  increases  severalfold 
in  crossing  the  shock.  Only  as  the  shock  nears  the  rigid  wall  do  the  zone 
masses  behind  the  shock  approach  their  initial  values;  in  the  von  Neumann- 
Richtmyer  scheme,  the  masses  do  not  change  at  all.  When  the  shock  gets 
close  to  the  rigid  wall,  the  zoning  in  the  unshocked  material  is  much  finer  than 
at  the  start  of  the  problem,  and  the  error  in  the  reflected  shock  pressure  is 
considerably  less  than  in  the  von  Neumann-Richtmyer  case. 

It  can  be  seen  that  there  are  some  very  simple  and  fundamental  problems 
which  have  not  been  solved  adequately  in  the  past  —  at  least  in  the  sense  that 
the  method  of  characteristics  affords  an  adequate  solution.  In  view  of  the 
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machine  time  spent  on  hydrodynamic  problems,  it  would  appear  that  a  good 
deal  more  investigation  of  the  difference  equations  themselves  is  in  order. 
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2  Technical  Library 

1  DCAS  (Tech  Ubrary),  AF  Unit  Post  Office,  Los  Angeles  45,  CmUf 

BSD,  AF  Unit  Post  Office,  Los  Angeles  45,  Calif 
1  (WDTV)  ATTN:  Col  Middlekauf 
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(Tech  Library) 

SSD  (Tech  library).  AF  Unit  Poet  Office,  Los  Angelee  45,  Calif 
ESD,  Hanscom  Fid,  Bedford,  Maes 
(Tech  Library) 

(ESAT) 

(ESF) 

AF  Mel  Dev  Cen  (MDGRT),  Holloman  AFB,  N  Mex 

European  Office  of  Research.  ATTN:  MaJ  Sullivan,  Brussels, 
Belgium 

KIR T LAND  AFB  ORGANIZATIONS 

AFSWC.  Kirtland  AFB,  N  Mex 

1  (SWNH) 

100  (SWOI) 

10  (SWRP) 

2  (SWRS) 

1  US  Naval  Weapons  Evaluation  FaciUty  (NWEF)  (Code  404), 
KirtUnd  AFB,  N  Mex 

OTHER  AIR  FORCE  AGENCIES 

Director,  USAF  Project  RAND,  via:  Air  Force  Liaison  Office, 
The  RAND  Corporation,  1700  Main  Street,  Santn  Monica,  Calif 

2  (RAND  Library) 

1  ATTN:  Dr.  O.  A.  Nance 

1  Dr.  A1  Latter 

1  Dr.  R.  L.  Bjork 

1  Dr.  Carl  Greifinger 

1  Mr.  Kenneth  Schwarts 

1  Lt  Col  Whitener 

ARMY  ACTIVITIES 

1  Chief  of  Research  and  Development,  Department  of  the  Army 

(Special  Weapons  and  Air  Defense  Division),  Wash  25,  DC 


No.  Cys 

2 

1 

1 

1 

1 

1 
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Director,  Ballistic  Research  Laboratories,  Aberdeen  Proving 
Ground,  Md 

1  Library 

1  ATTN:  Mr.  Frank  J.  Allen 

1  Mr.  E.  D.  Baicy 

1  ARGMA  (Tech  Library),  Huntsville,  Ala 

NAVY  ACTIVITIES 

I  Chief  of  Naval  Operations,  Department  of  the  Navy,  ATTN:  OP-36, 

Wash  25,  DC 

Commanding  Officer,  Naval  Research  Laboratory,  Wash  25,  DC 
1  (Tech  Library) 

1  ATTN:  Mr.  Walter  Atkins 

1  Mr.  Carroll  Porter 

1  Commanding  Officer,  Naval  Radiological  Defense  Laboratory 

(Technical  Info  Div),  San  Francisco  24,  Calif 


OTHER  DOD  ACTIVITIES 

Chief,  Defense  Atomic  Support  Agency,  Wash  25,  DC 

1  (Document  Library) 

1  ATTN:  Lt  Col  Singer 

1  Maj  McCormac 

Commander,  Field  Command,  Defense  Atomic  Support  Agency, 
Sandia  Base,  N  Mex 

1  (FCAG3,  Special  Weapons  Publication  Distribution) 

1  ATTN:  Lt  Col  Morgan 

Director,  Advanced  Research  Projects  Agency,  Department  of 
Defense,  The  Pentagon,  Wash  25,  DC 

1  ATTN;  Dr.  Charles  Cook 

1  Lt  Col  Roy  Weidler 

1  Director,  Defense  Research  &  Engineering,  ATTN:  Col  Gilbert, 

The  Pentagon,  Wash  25,  DC 

10  ASTIA  (TIPDR),  Arlington  Hall  Sta,  Arlington  12,  Va 
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AEG  ACTIVITIES 

US  Atomic  Energy  Commiseion  (Technicnl  Report*  library,  Mr*. 

J.  0"Leary  for  DMA).  Wa*h  25.  DC 

Preaident,  Sandia  Corporation,  Sandia  Ba**,  N  Mex 
(Document  Control  Division) 

(ATTN:  Mr.  C.  O.  Lundergan) 

Director,  University  of  California  Lawrence  Radiation  Laboratory, 

Technical  Information  Division,  P.  O.  Box  808,  Livarmore,  CalU 

ATTN:  Mr.  Clovis  Craig 
Dr.  John  G.  Trulio 
Dr.  Mark  Wilkins 
Mr.  Robert  Banaugh 
Mr.  Joseph  Fleck 
Mr.  James  Quong 

Director,  Los  Alamos  Scientific  Laboratory,  P.  O.  Box  1663, 

Lo*  Alamos,  N  Mex 

(Helen  Redman,  Report  Library) 

(J-15,  Dr.  Arthur  Cox) 

(T-5,  Dr.  George  White) 

(J-10,  Dr.  Herman  Hoerlin) 


OTHER 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


Edgerton,  Germeshausen,  &  Grier,  ATTN:  Dr.  Marion  Schuler, 

150  Brookline  Ave. ,  Boston.  Mass 

General  Dynamics/Convair,  Physics  Division,  ATTN:  Dr.  C.  G. 
Davis,  Box  1950,  San  Diego  12,  Calif 

Cornell  Aeronautical  Lab,  4455  Genessee  Street,  Buffalo  21,  NY 
ATTN:  Dr.  Walter  Gibson 
Dr.  J.  Gordon  Hall 

Institute,  of  Aerospace  Sciences,  Inc.,  2  E.  64th  St. ,  New  York,  NY 
OTS,  Department  of  Commerce,  Wash  25,  DC 
Stanford  Research  Institute,  Menlo  Park,  Calif 
ATTN:  Dr.  C.  M.  Ablow 

Dr.  George  Abrahamson 
Mr.  John  O.  Erkman 
Mr.  G.  R.  Fowles 
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1 

1 
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1 


1 

1 

1 
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Aeronutronic,  A  Division  of  Ford,  Ford  Road.  Newport  Beach,. 
Calif 

ATTN:  Dr.  Raymond  Grandey 
Dr.  Donald  Sachs 

Mr.  John  W.  Lipford  . 

The  Boeing  Company,  Aerospace  Division,  Box  3707,  Seattle  24, 
Wash 

ATTN:  Dr.  Donald  Hicks 
Dr.  Donald  Keller 
Mr.  Daniel  Young 

Massachusetts  Institute  of  Technology,  Aeroelastic  and  Structures 
Laboratory,  ATTN:  Dr.  Walter  Herrmann,  Cambridge  39,  Mass 

Aerojet  General  Corporation,  1711  Woodruff  Ave. ,  Downey,  Cmlif 

ATTN:  Mr.  Nigel  Thomas 

Mr.  Mark  Wagner 

Dr.  Louis  Zemow 

General  Dynamics  /  General  Atomic,  Box  608  San  Diego  12,  Calif 
ATTN:  Dr.  Burton  Freeman 
Dr.  Charles  Loomis 

RCA,  Advanced  Military  Systems  Group,  David  Samoff  Research 
Center,  Princeton,  NJ 

ATTN:  Mr.  Fred  Herzfeld 

Dr.  John  Jarem 

Geophysics  Corporation  of  America,  ATTN:  Dr.  Harry  £.  Stubbs, 
Bedford,  Mass 

AVCO-RAD,  ATTN:  Dr.  Jerrold  Yos,  201  Lowell  Street, 
Wilmington,  Mass 

Lockheed  Missiles  and  Space  Division,  3251  Hanover  Street, 

Palo  Alto,  CaUf 

ATTN:  Dr.  R.  E.  Meyerott 

Dr.  R.  K.  M.  Landshoff 

Dr.  Dan  Holland 

Republic  Aviation  Corporation,  Plasma  Propulsion  Laboratory, 
Farmingdale,  NY 

ATTN:  Mr.  Milton  Halem 

Dr.  William  Mcllroy 


191 


TN-61-29 


No.  G.Yi 


DISTRIBUTION  (con't) 


1  McAllister  and  Associates,  ATTN:  Dr.  Louis  E.  Bothell,  203 

Truman,  N.  E. ,  Albuquerque,  N  Mex 

1  Allied  Research  Associates,  ATTN:  Mr.  Sheldon  Kahalas,  43 

Leon  Street,  Boston,  Mass 

1  General  Electric  Company,  Aeroscience  Laboratory,  ATTN:  Dr. 

F.  A.  Lucy,  3750  "D”  Street,  Philadelphia  24,  Pa 

Stanford  Research  Institute,  Menlo  Park,  Calif 

1  ATTN:  Dr.  L.  Evan  Bailey 

1  Dr.  Ernest  G.  Chilton 

1  American  Science  and  Engineering,  ATTN:  Mr.  Gilbert  Fryklund, 

79  Broadway,  Cambridge,  Mass 

Space  Technology  Labs,  Box  95001,  Los  Atigeles  45,  Calif 
1  ATTN:  Dr.  H.  Leon 

1  Dr.  J.  Maxey 

1  California  Institute  of  Technology,  ATTN:  Dr.  Max  Williams, 

Pasadena,  Calif 

1  Battelle  Memorial  Institute,  ATTN:  Mr.  B.  J.  Brand,  505  Kix^f 

Avenue,  Columbus,  Ohio 

1  Technical  Operations,  Inc. ,  ATTN:  Dr.  Paul  I.  Richards, 

Burlington,  Mass 

1  Official  Record  Copy  (SWRPA,  Lt  Rich) 
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